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AERODYNAMIC PRELIMINARY ANALYSIS SYSTEM II 


PART I THEORY 


By E. Bonner, W. Clever, K. Dunn 


North American Aircraft Operations, Rockwell International 


SUMMARY 


An aerodynamic analysis system based on potential theory at subsonic/ 
supersonic speeds and impact type finite element solutions at hypersonic 
conditions is described. Three-dimensional configurations having multiple 
non-planar surfaces of arbitrary planform and bodies of non-circular contour 
may be analyzed. Static, rotary, and control longitudinal and lateral- 
directional characteristics may be generated. 


The analysis has been implemented on a time sharing system in 
conjunction with an input tablet digitizer and an interactive graphics 
input/output display and editing terminal to maximize its responsiveness to 
the preliminary analysis problem. Computation times on an IBM 3081 are 
typically less than one minute of CPU/Mach number at subsonic, supersonic or 
hypersonic speeds. Computation times on PRIME 850 or a VAX 11/785 are about 
fifteen times longer than on the IBM. The program provides an efficient 
analysis for systematically performing various aerodynamic configuration 
tradeoff and evaluation studies. 
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INTRODUCTION 


Aerodynamic numerical analysis has developed to a point where 
evaluation of complete aircraft configurations by a single program is 
possible. Programs designed for this purpose in fact currently exist, but 
are limited in scope and abound with subtleties requiring the user to be 
highly experienced. Many of the difficulties are attributable to the 
numerical sensitivity of the associated solution. In preliminary design 
stages, some degree of approximation is acceptable in the interest of modest 
turn-around time, reduced computational costs, simplification of input, and 
stability and generality of results. The importance of short elapsed time 
stems from the necessity to systematically survey a large number of 
candidate advanced configurations or major component geometric parameters in 
a timely manner. Modest computational cost allows a greater number of 
configurations and/or conditions to be economically investigated. 

One approach in this spirit is to employ panel approximations which 
reduce the number of simultaneous equations required to satisfy flow 
boundary conditions. Surface chord plane formulations, locally two 
dimensional crossflow body solutions and non - in t e r f er ing panel 
simplifications are examples of approximations which can be used for this 
purpose. An alternative approach is to use surface chord plane formulations 
again for thin surfaces which can carry lift and surface panels for thick 
body type regions . 

Finite element analysis when combined with realistic assessment of 
limitations and estimated viscous characteristics provides a valuable tool 
for analyzing general aircraft configurations and aerodynamic interactions 
at modest attitudes for subsonic/supersonic speeds and evaluation of 
compressible non-linearities at high Mach numbers. 
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LIST OF SYMBOLS 


A 

A. . 


~AVG 


c 0 , c , c 

2 m n 


n 


pNET 


Projected oblique cross -section area 

Influence coefficient. Normalwash at control point i due 
to vortex panel j of unit strength 

Area of quadrilateral panel i 

Reference span 
Local chord 
Reference chord 
Average chord 

Section drag coefficient 

Drag coefficient 

Flat plate average skin friction coefficient 
Boundary condition for control point i 
Section lift coefficient 

Rolling, pitching and yawing moment coefficients 
Lift coefficient 

Section normal force coefficient 
Pressure coefficient (P-P )/q 

Net pressure coefficient (P^-P )/q and vortex panel strength 
Leading edge suction coefficient 
Leading edge thrust coefficient 


C ,C ,C Axial, side, normal force coefficient 
x y z 
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F’x’Fy’F’z Force components 


g(x) Axisymmetric outer solution to potential equation 

h Radius of curvature of cross-sectional boundary 


i, j ,k 
K 

K 

s 


i 

i(i,n) 


L 

L/d 

M 


M ,M ,M 
x y z 


Unit vectors in x,y,z direction respectively 

Drag due to lift factor or skin friction thickness correction 
factor 

Equivalent distributed sand grain height or attainable suction 
fraction 

Effective length 

Length of segment i, i+1 of contour 

Equivalent body length or geometric length 
Body fineness ratio 
Mach number 
Moment components 


n Unit normal 

p,q,r Rolling, pitching and yawing velocity about x, y and z 


P.q>r 



Nondimens ional angular velocities pb/2U, qc/2U and rb/2U 
Static pressure 
Prandtl number 

2 

Free stream dynamic pressure, l/(2pU ) 

Recovery factor 

Unit Reynolds number or radius of curvature 
Reynolds number based on [ ] 


R 


Gas constant 



^ to 


s 

S 

S REF 

T 


t/c 

u,v,w 

u,v 

V. 

J 

W 


x, y,z 
x f r, tf 
Z 
a 
a. 


r 

*? 

8 . . 
ij 


8 . 
J 



5i//£x 


A 


Segment arc length 

Body cross-sectional area or surface area 
Reference area 

o 

Static temperature , R, or tangent of quadrilateral panel 
leading edge sweep 

Airfoil thickness ratio 

x,y,z nondimensional perturbation velocity components 
Freestream velocity 
Jet velocity 

Complex potential function 

Body axis Cartesian coordinate system 

Body axis Cylindrical coordinate system 

Complex number y+iz 

Angle of attack 

Local angle of attack at surface control point i 
Angle of sideslip or [|l-M | ] 

Vorticity strength per unit length or ratio of specific heats 

Horseshoe vortex strength in Trefftz plane 

Deflection or impact angle 

Lateral surface coordinate 

Kroneker delta: 0 i^j 

1 i-j 

Jet deflection angle relative to trailing edge 
Total jet deflection angle 
Body slope 

Dihedral angle of quadrilateral panel or boundary layer momentum 
thickness 

Sweep angle 
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T 

</> 

$ 

Q 

Subscript 

c 

CG 

e 

F 

£ 

LE 
r 
t 
T 

TRAN 
u 
v 
w 


Absolute viscosity 

Kinematic viscosity, p/p 

Density 

Source density 

Side edge rotation factor 

Perturbation velocity potential 

Total velocity potential 

See figure 4 

Leading edge rotation factor 
camber 

center of gravity 

edge conditions 

friction 

lower surface 

leading edge 

recovery 

thickness 

tip 

transition point 
upper surface 
vortex 
wave 

freestream condition 


Superscripts 

first derivative or quantity based on effective origin 
M second derivative 

* Eckert reference temperature condition 

■+ vector quantity 
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SUBSONIC/SUPERSONIC 


The arbitrary configurations which may be treated by the analysis are 
simulated by a distribution of source and vortex singularities. Each of 
these singularities satisfies the linearized small perturbation potential 
equation of motion 


2 

fi <t> 


+ d> + dt — 0 
xx yy zz 


The singularity strengths are obtained by satisfying the condition that 
the flow is tangent to the local surface : 


d$/dn = 0 

All of the resulting velocities and pressures throughout the flow may be 
obtained when the singularity strengths are known. A configuration is 
composed of bodies, interference shells and aerodynamic surfaces (wings, 
canards, tails etc.). There are two alternative types of singularities used 
to represent the configuration. Figure 1 shows the first type, which can be 
used at all Mach numbers, and figure 2 shows an alternative method, which 
can be used only at subsonic Mach numbers. 


wing and vertical tail 



Figure 1 A. Singularities Used to Simulate a Configuration. 
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In the first method, the first step in the solution procedure consists 
of obtaining the strengths of the singularities simulating the fuselage and 
nacelles, from an isolated body solution. The present analysis uses 
slender-body theory to predict the surface and near field properties. The 
solution is composed of a compressible axisymmetric component for a body of 
revolution of the same cross-sectional area and an incompressible crossflow 
component, <f> , satisfying the local three dimensional boundary conditions in 
the (y,z) plane. The crossflow is a solution of Laplace's equation 


9 


yy 


+ 


zz 


0 


A two-dimensional surface source distribution formulation is used to obtain 
this solution. When the body singularity strengths are determined, the 
perturbation velocities which they induce on the aerodynamic surfaces , or 
other regions of the field, are evaluated. 


The assumptions of thin airfoil theory allow the effects of thickness 
and lift on aerodynamic surfaces to be considered independently. Therefore, 
the effects of the aerodynamic surfaces can be simulated by source and 
vortex singularities accounting for the effects of thickness and lift, 
respectively. The source and vortex distributions used in this program are 
in the form of quadrilateral panels having a constant source or vortex 
strength. The vortex panels have a system of trailing vorticies extending 
undeflected to downstream infinity. The use of a chordwise linearly varying 
source panel is provided as an a option to eliminate singularities 
associated with sonic panel edges at supersonic Mach numbers. The panels 
are planar, that is they have no incidence to the free stream (although 
dihedral may be included) , since thin airfoil theory allows the transfer of 
the singularities and boundary conditions to the plane of the mean chord. 
These boundary conditions are satisfied at a single control point on each 
panel. For thickness, the control point is located at the panel centroid 
while the effects of twist, camber, and angle of attack are satisfied at the 
spanwise centroid of each vortex panel and at 87.5 percent of its chord. 


A cylindrical, non-circular, interference shell, composed entirely of 
vortex panels, is used to account for the interference effects of the 
aerodynamic surfaces on the fuselage and nacelles. The boundary conditions 
on an interference shell are such that the velocity normal to the shell 
induced by all singularities, except those of the body which it surrounds, 
is zero. The boundary conditions are satisfied at the usual control points 
for vortex panels. 


The second alternative method uses constant doublet panels and constant 
source panels to represent the body surface. These panels can be of an 
arbitrary quadrilateral shape and may be inclined to the direction of flow. 
The aerodynamic surfaces are represented by the same type of chord plane 
source and vortex panels as were used in the first method. 
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Alternative method for subsonic flow only 


wing and vertical tail 
- chord plane source and vortex panels - 



Figure 1 B. Singularities Used to Simulate a Configuration (M qq < 1), 

This second method can be used at subsonic Mach numbers only . At 
supersonic Mach numbers, the doublet panels, which are equivalent to 
quadrilateral vorticies, produce infinite perturbation velocities in certain 
regions of the flow, and thus cannot be used. The body source and doublet 
strengths are chosen to satisfy both an arbitrary normal velocity boundary 
condition on the body, 


d$/dn - V n 

and to have zero perturbation potential in the entire region interior to the 
body surface . 


$ - 0 


The following sections define the details of the solution procedure. 
Included are discussions of the isolated body analysis, surface finite 
element analysis considering edge effects, and evaluation of aerodynamic 
characteristics including drag. References are cited for the reader 
interested in further pursuing a particular point. 


8 



SLENDER BODY SOLUTION 


1 2 

According to slender body theory 9 the flow disturbance near a 
sufficiently regular three-dimensional body may be represented by a 
perturbation potential of the form 

<t> - <p(y,z;x) + g(x) (1) 

<p(y,z;x) is a solution of the 2-D Laplace equation in the (y,z) cross flow 
plane satisfying the following boundary conditions 

V<p“jv + kw - 0 

d$/dn - 0, on C(x) (2) 

C(x) and n are defined in figure 2. A general solution for <f> may be written 
as the real part of a complex potential function W(Z) with Z - y+iz . 


<p « R g W - R & [ A 0 (x) In Z + 2 A n < x ) Z " n 5 

n=l 


A useful alternative representation of <f> and W is obtainable with the aid of 

3 

Green' s theorem . 


- R W « -2 R i a(D ln(Z-f)ds (3) 

C(x) 

where a(f) is a "source" density for values of f - y + iz c , (y c » z c ) being 
coordinates of a point on the contour C(x). 

The function g(x) obtained by matching <f> of equation (1) which is valid 
in the neighborhood of the body with an appropriate "outer" solution. g(x) 
is then found to depend explicity on the Mach number M and longitudinal 
variation of cross-sectional areas S(x) 

g(x) - l/(2ir)[S'(x)ln(O.50)-l/2f S" (t)ln(x- t)dt -hl/2 f S" ( t) ln( t-x)dt 

J 0 J x 

-1/2 S' (0) In x -1/2 S' (l)ln(l-x)] M < 1 

(4) 

g(x) - 1/(2 tt)[ S ' (x) ln(0 . 5/3) - f S" ( t)ln(x- t)dt ] M > 1 

J 0 

The body axis perturbation velocities are obtained by di f f er intiat ion 
of equation (1) 
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y S (ds)=5 









V — 


u “ - <P X + g' (x) 


w = <j> 


At supersonic speeds, zone of influence considerations require that u - 
w = 0 for x - p r < 0. 

Solution of the preceding equations is based on an extension of the 
method of reference 3. 


CROSS FLOW COMPONENT 

The reduction of computations to a numerical procedure utilizes the 
integral representation of <p given in equation (3) by discretization of the 
cross-sectional boundary into a large number of short linear segments 
(figure 3) over each of which the source density a is assumed constant at a 
value determined by boundary conditions. 

Computation of a(i,n) over the segment i, i+1 proceeds by applying the 

boundary condition equation (2) at each segment of C . If V<p - q - jv + icw 

represents the velocity vector, the corresponding complex velocity in the 
cross flow plane is obtained by differentiation of W in equation (3) with 
respect to 2: 


v - iw = -2j) a(f)/( Z-f) ds (5) 

The contribution by the sources located on segment i, i+1 to the velocity at 
^ j > n ^ rst evaluated. Noting that i, i+1 makes an angle 0(i,n) with 

respect to the horizontal axis, we have 


, i0(i,n) 

df * ds e v ' 

and the contribution ot the integral in equation (5) may be written: 

f i+l,n 

A[v(j,n) - iw(j ,n) ] - -2a(i ,n)e' 10 (l ,n) f [Z. 

J J 


l , n 


11 




After integration of the last term and summation over all contributing 
segments, the result may be written 

v(j ,n)-iw(j ,n) — 2^ a(i,n)e ^ ^ ,n ^ln[R( i+1 , j ,n)/R(i , j ,n) ]+iS (i , j ,n)^ (6) 

i 

in which referring to figure 4, the quantities R(i,j,n) and 6(i,j,n) are 
defined by the relationships 


R(i » j ,n)e 


iV>(i, j ,n) 


Z. - 
J 


h,n 


5(i, j , n) - VK i , j ,n) - V>(i.j .n) 

To insure uniqueness of the complex velocity, care must be exercised in 

assigning values to the angles ^(i,j,n) and^(i,j,n). Referring to figure 
4, these are measured counter-clockwise from the positive y-axis so that 

when facing P. to P . , a point P just to the left of i , i+1 shall 

define an angle = 0(i,n). As P traverses a path around 

J > 11 

P to a point just to the right of i,i+l, ^(i,j,n) increases from 0(i,n) 
i ,n 

to 0(i,n) + 27 t. The same holds true for V>(i,j,n) as P. traverses a path 
around 

P In consequence of these definitions S(i,j,n) becomes -n when 

i + l,n 

approaching i,i+l from the right and n when approaching from the left. This 
discontinuity reflects that exhibited by the stream function upon traversing 
any closed path which encloses a distribution of finite sources . 

From the boundary condition equation (2) , we have 

-(dip/d n) . = v(j ,n)sin0 (j ,n) - w(j ,n)cos0 ( j ,n) 

J 

After substitution of v and w from equation (6), this last expression 
becomes 

-(dcp/dn). = ) a(j , i)a(i.n) (7) 

J j n L 


where 


a(j , i) 


2-|sin[ 0 ( j ,n) - 0<i,n)] ln[R(i+l,j ,n)/R(i, j ,n) ] 
+ 5 (i , j ,n) cos [ 0 ( j ,n) - 0(i,n)]\ 
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j ,n) 

Figure 4. Details of Variables Pertaining to Segment i,i+l of Boundary C . 

n 
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The surface normal perturbation velocity - (dtp/dn) ^ ^ may be written in 

terms of the body slope (dv/dx) . , the angles of attack a, and sideslip 0 

J > n 

and the angular velocities p,q,r as 


- (dcp/d n) 

J .n 


( dv/dx ) + [a + q(x-x )/U + py/U] cos0(j,n) 

J > 11 C & 


+ [0 - r(x-x )/U + p(z-z )/U] sin0(j,n) 
C S c 6 


Satisfying equation 7 at each of the points 
yields a set of equations for a(i,n). 


P. 

J 


on a given contour boundary 


AXISYMMETRIC COMPONENT 

Differentiation of g(x) must be carried out with due concern for the 
nature of the improper integrals appearing in equation (4) . The result is 

g'(x n ) - l/(4w) js"(x n )ln[0.25(l-M 2 )] + I n (x n ) - J n (x n ) 

-S'(0)/x + S' (l)/(l-x ) - S"(0)ln x - S" (l)ln(l-x )l M < 1 

n n n n n J 

g'(x ) - l/(2ir) -j 1/2 S" (x )ln[0.25(M 2 -l)] - J (x ) - S"(0)ln x \ M > 1 
n ^ n n n nj 

where 

l n-1 

I ri ( x _ ) - S ln(x -t) S'"(t)dt = ) [S- - S” ] ln(x -x ) 

n n ^ n (_ m+1 m J m n 

n m=n 

x 

n 

J (x ) - f ln(x -t) S' " (t)dt 
n n J 0 n ' v ' 

x - (x -+ x )/2 
m m+1 m 


n-1 

- > [S- - S'* ] ln(x -x ) 

l_ m+1 m n m 

m=0 


To compute the second derivatives of the equivalent body cross-sectional 

area required for g' (x) , the first derivatives at x are found by finite 

m J 


differences between x and x Second derivatives S M (x ) 

m m+1 x m 

x )/2 are then found by finite differences between S' at 
Finally S"(x m ) is determined by linear interpolation of S ” 


at x -(x - + 

m m+1 

x and x ^ . 
m m+1 


( x ) between 
m 


x 

m 


and x 


m+1 * 
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PERTURBATION VELOCITIES 


t f ^ ie ax i a ^ velocity u depends on dcp/dx and the axisymmetric solution 
g'(x). dcp/dx is obtained by differentiation of the integral in equation (3) 
obtain an exact expression which is then approximated by evaluating 
the result over the segmented boundary. 

The derivation of dcp/dx must take into account the fact that the path of 
integration in equation (3) is a function of x. Referring to figure 2 
increments of a dependent variable taken along C(x) are denoted by d( ) and 
increments taken normal to C are denoted by 6( ). Differentiation of 
equation (3) then yields 

dcp/dx - -2 R e [ | ( Sa/Sx ) ln(Z-f)ds - | a(r)/(Z-f) (fif/«x)ds 

+ | a(0 ln(Z-f) (8(ds)/8x ) ] (8) 

From figure 2 

6(ds) - 8 vd$ - Si j ds/h(f) (9) 

where h(f) is the radius of curvature of C(x) at f. In addition, we have 
from figure 2 

6$/6x - 8u/8x e^^ 0.5 jt) (10) 

To evaluate 8o/8x we note, 


6a/6x - lim [a(i,n+l) - a(i,n)]/$x (11) 

8x-> o 

Introducing equations (9), (10), and (11) into equation (8), 

dcp/dx -= -2 R^ |[(5o/5x) 0 + a / h 8v/8x ] ln(Z-f)ds + i j)[o(5i//Sx) ] df/(Z-$-)j- 

Again, assuming that quantities in the brackets of the integrands are 
constant over i,i+l, 

W/3X) j >n " 2 } {t(W«x) 0 + o/h(6u/Sx)] i n A<p(i, j ,n)/a(i,n) 
i 

o(i,n) (8v/8x) i ^ 5(i , j ,n)J 
| R(i+l,j,n).u(i,n) In R(i+l,j,n) 

- R(i, j ,n) »u(i,n) In R(i,j,n) 

- R(i,j ,n)*n(i,n) 6(i,j,n) + i(i,n) j- 

The radius of curvature h(i,n) and the derivatives 8a/8x, 
approximated at the mid-points of the segments i,i+l as follows 

16 


where 

A<p(i, j ,n)/cr(i,n) = 


8v/8x are 



8a/8x - the derivative at the mid-point x^ of the interval 

x , x , is set equal to the divided difference between 
n n+1 n 

a(i,n) and a(i,n+l). Linear interpolation between these 


1 A 


derivatives then yields 8o/8x at x 


n 


8u/8x - referring to figure 5, the displacement 8t} is determined 
by linear interpolation between ^ and n * 

^/( X n+r X n) then represents 8v/8x at x^. Linear 
interpolation between the stations x'^ then yields 8v/8x 


at x . 
n 

6 at ^i n determined by interpolation between values 

of 0(i,n) at P. . The curvature 1/h at P. is then set 
i,n ' i , n 

equal to the divided difference between 6 at P, 1 and 6 

^ i+l, n 

at P. 

i,n. 


The lateral and vertical perturbation velocities, v and w , are obtained 

from 

v - iw= - 2 j)a(f)/(Z - D ds 

Integration over the boundary with constant segment source density yields: 
v(j,n) - iw (j ,n) - 2^ o(i ,n)e id ^ i ,T ^jln[R(i+l , j ,n)/R(i , j ,n) ] - i5(i,j,n)j- 


Thus 

V = <f> 


^ a(i,n)| 


ln[R(i+l, j ,n)/R(i, j ,n) ]cos 8 ( i,n) - fi(i,j,n) sin 0(i 


i,n)l 


w - <j> z - 2 ^ {7(i,n)^ln[R(i+l, j ,n)/R(i, j ,n) ]sin 0(i,n) - 6(i,j,n) cos 0(i,n)j- 
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SURFACE SOLUTION 


CHORD PLANE SOURCE AND VORTEX PANELS 


The wing, canard, vertical and horizontal tail are simulated by a 
system of swept tapered chord plane source and vortex panels with two edges 
parallel to the free stream. The coordinates of the panel corners are 
specified with respect to an (x,y,z) system having its x-axis in the free 
stream direction and its z-axis in the lift direction. The panel influence 
equations are written in terms of a coordinate system having a z-axis normal 
to the panel and an x-axis along one of the two parallel edges. A 
coordinate transformation is necessary to obtain the coordinates in the 
panel reference system. If the plane of the panel is inclined at an angle 
^p res P ect to the y,z plane, a transformation into the panel coordinate 


system (x ,y ,z ) 
P ; P P 


is accomplished as shown in figure 6. 


x - x 
P 


y - y 
P 

cos 

e 

p 

+ z 

sin 

8 

P 

z - - y 
P 

sin 

8 

P 

+ z 

cos 

8 

P 




control 

point 

panel 


w 

c 



u 

c 


V 

c 


w 

c 


u 

p 


“ V cos 
p 

O 

l 

V 

+ 

w sin 
P 

<V 

*c> 

= -v sin 
P 

l 

o 

V 

+ 

w cos 
p 

<V 

V 


Figure 6. Coordinate Transformation in Panel Reference System. 
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A transformation of the (u ,v ,w ) velocities into the coordinate system of 

P P P y 

the panel on which the control point is located (u ,v ,w ) results in the 

c c c 

axial, binormal and normal velocities induced on the panel. 

For the image of the influencing panel, the signs of y .6 and v are 

c c c 

changed while using the same calculation procedure. 


QUADRILATERAL SOURCE AND DOUBLET PANELS (PANELED BODIES) 


For subsonic Mach numbers the body may be represented by a system of 
planar quadrilateral constant source and constant doublet panels. Since 
four points arbitrarly selected on a surface may not lie in the same plane, 
a mean surface through the four points is selected to represent the panel. 

Let (x^,y^,z^) represent the four points on the body surface, 


and ^i ,r *i’^i^ represent the four points on the mean surface. 

(x,,y 2 .2 2 ) 


> ^2 > T2 ) 



(^3.^ 3 .r 3 ) 

(X3 , y 3 , 2 3 ) 


This mean surface is chosen in the following manner. 

1. The direction of the panel normal is found by taking the cross product 
of the vectors representing the diagonals. 

_ d 3i x d 42 d 31 - ( x 3 - Xl , y 3 - yi , Zj-Zj ) 

11 “ ” 

I d 3ll I d 4 2 I d 42 - ( X 4‘ X 2 > > Z 4‘ Z 2 ) 
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2. The out of plane distance, 6 t is calculated using vectors determined by 
pairs of points. 


( s 12 + s 34 ) 


n 


( x 1 -x 2 , y x -y 2 , z x -z 2 ) 

( x 3 -x 4 , y 3 -y 4 , z 3 -z 4 ) 

The coordinates of the mean surface are calculated by adding or 


*1 2 


>3 4 


subracting 8 

-► 

n 

from 

each 

of 

the corner 

points 

i 

. e . 

( *i. 

h • 

fi> 

- 

( 

x i . 

■ yi 

. Zj) 

- 

S 

( n i » 

n 2 , 

n 3 ) 

( e 2 . 

I 2 . 

5*2 ) 

- 

( 

x 2 , 

■ y2 

. z 2 ) 

+ 

S 

( , 

n 2 , 

n 3 ) 

( £ 3 . 

1 3 , 

5*3 ) 

- 

( 

*3 - 

y 3 

. z 3 ) 

- 

6 

( , 

n 2 , 

n 3 ) 

( £ 4 . 

14 • 

5*4) 

- 

< 

x 4 , 

y 4 

■ Z 4 ) 

+ 

S 

( , 

n 2 » 

n 3 ) 


The normal computed for these four points is the same as the normal for the 
original body points, since the diagonal vectors are the same. If a vector 

determined by the line segment joining any two of the four points is normal 
-+ 

to n, then the four points must lie in the same plane. This is easily shown 
to be true. From the above definitions, 


( °12 + ^34) * n 


1 2 


'34 


( £ 1 ~ £2 > r ll~ r l2 ’ f 1 * f 2 ) 
( £ 3 -£ 4 . 13-14 , f 3 - f 4 ) 


but 


r 34 ” < £3^4 - 13-14 > C3-C4 ) 

- [ <£ 3 -£l) - «4-<2> + (£l*£ 2 ) 

. (la'll) - ( 14 - 12 ) + (h- 12 ) 

' ( 5" 3 ~ 5* 1 ) ' (r 4 -f 2 ) (fi-fa) 


or 


( a i 2 + a 3 4 ) • n - 2 o 12 * n = 0 

similarly 


d 3 1 - d 42 + °1- 


°\ 2 ' n 
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and all four points must lie in the same plane. 
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BOUNDARY CONDITIONS 


Vortex Panels 


At the vortex panel control points the resultant velocity along the normal 

at a panel control point must be zero. Using a local coordinate system, 

with perturbation velocities (u ,v ,w ) at the control points, 

c c c 



U*n =* U [ e + u e + (v + a sin 6 ) e, + (w + a cos 6 ) e 1 *n 
x c x c cb c c n 


= U [ (1+u )e *n + (w + a cos $ ) 1 = 0 
c X c c 


with n - [e - (dz /dx^e ] 
n c' 'i x J 

For small perturbations (1+u )e *n ^ - (dz /dx) . and e*n ~ 1 

c x c i n 


Therefore w * (dz /dx) . - a cos $ 

c 7 'i c 


Body Panels 

The boundary condition on body panels will involve the normal component of 
velocity. If we set the normal component equal to zero, we have the usual 
flow tangency boundary condition. Nonzero normal components can be used for 
jets or inlets. Given the boundary condition on the surface and the field 
at infinity, the solution for the external flow is unique. It can be 
satisfied by an infinite number of combinations of source and doublet 
distributions on the surface. However each combination will result in a 
different field inside the body surface. Specifying an internal boundary 
condition will make the source and doublet distribution unique, and can have 
a powerful effect on the numerical behavior of a solution involving a finite 
number of elements. The internal boundary condition which we have chosen, 
with these numerical considerations in mind, has zero perturbation potential 
on the internal boundary, and therefore due to the nature of the governing 
equation, zero perturbation potential inside. Below, we will show that by 
first correctly choosing the surface source distribution, we can also 
satisfy the external normal velocity boundary condition by satisfying the 
internal surface boundary condition on <f > . 

Consider a closed region determined by the surface S. Let the 
surface have a distribution of sources and doublets with local strength o 
and fi. The surface, S, will divide the interior and exterior regions. 
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exterior e 




the subscript denoting properties exterior to S 
the subscript denoting properties interior to S 


the external normal to S » ( n n n) 

x y, z 


V *= perturbation velocities due to a and \x . 


the prescribed normal velocity on S (exterior) 


We can set the value of the surface source strengths to any value and still 
satisfy the external boundary condition. 

A _> o 

We will set, a = - u • n + U 

oo e n 


and adjust the value of and any other singulartity strengths, such that 
everywhere on the interior surface of S , 

4 >. = o 

i 

Then in the entire region interior to S, 

4 >. - o 

i 




and 



U. - V <f> + U - U 

1 i CO 00 

Since the p gives a continuous normal velocity across S, using Appendix C, 


u • n + u. • n # = u-n 
e e l i e e 


- U • n + U 
oo e n 


( U + u ) • n 
oo e e 


as required on S . 


Therefore the normal velocity boundary condition can be satisfied by 
substituting a boundary condition for <f> on the internal boundary surface . 
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PANEL INFLUENCE COEFFICIENTS 


Each of the panel types induces a perturbation potential everywhere in 
space. If panel j has unit strength, we can say it will induce the 
following velocities and velocity potential at the control point of panel i. 


(A U . , 

A V . , 

A W ., 

A?.) 

for vortex panels 

(D?. , 

D ij' 

D. . , 
iJ 

of.) 

1J 

for body doublet panels 

(S?. , 

»Ir 

S w ., 

1J 

sf.) 

1J 

for body source panels 

(T? . , 


T W ., 

T^ . ) 

for thickness source panels 


Therefore, assuming there are ntv vortex panels and ntb body panels, and 
C , m-, ° ■ . T ■ . are the P anel singularity strengths, the following set of 
Pj j j J 

panel influence coefficients can be written: 


ntv 


ntb 

ntb 


ntv 



u. = r c 

j-i J 

+ 

Y dY. h. + 
Z_ 3-J J 

r* 

a . + 
J 

r * u - 

Z_ 

r . + 

J 

U 0 

i 


j“l 

j=i 


j-1 


ntv 


ntb 

ntb 


ntv 



v, - Y~ aY. C 
i Z_ Pi 

+ 

r °T- /*• + 

/ ij j 


a. + 
J 

xx 

r . + 
J 

v o 

i 


j-i 

j-1 


j-i 


ntv 


ntb 

ntb 


ntv 



w. - Y~ aV. C 
i /_ Pi 

j“l J 

+ 

r ^ + 

r s -. 

L- 

a . + 

J 


r + 
P 

w o 

i 


j-i 

j-i 


j-i 


ntv 


ntb 

ntb 


ntv 



<j>. = ^ A?. C 

*1 L_ 1J P< 

j-1 J 

+ 

r of. m + 

Z_ P 

Y st. 

Z_ 

a + 
P 

Y t ^- 

L _ ij 

T + 

P 

<t>0 

i 





j-i 


where (u 0 , v 0 , w 0 , 

K 

) refer to the perturbations induced by any other 

i i i 


i 






body or source singularities, e.g. slender bodies. 





24 



PANEL SINGULARITY STRENGTHS 


Source (Thickness) Panels 

The source singularity strengths for thickness panels may be found 
directly by equating each source panel strength to the slope of the 
thickness distribution at its control point. For panel i 

'i - «V dx, i 

w ^ ei ^ e refers to the shape of the thickness distribution. 


Body Source Panels 


The source singularity strengths for body source panels are set to give 
the correct normal velocity boundary condition when the internal 
perturbation potential is zero. For panel i we set 


o. = - U • n. + U 

1 00 l n. 

l 


<7 . 
1 


2 2 2 2 ~1 
3 n + n + n 

X i z i J 


1/2 


where n. = ( n 

i x. 


n 


n z ) is the outward normal of panel i, and U is 


the normal velocity boundary condition for panel i. 


n. 

l 


Vortex and Body Doublet Panels 


The determination of the vortex and doublet panel singularity strengths 
is the final step in the solution procedure. They are obtained by solving a 
set of simultaneous equations utilizing the panel influence equations to 
relate the singularity strengths to the boundary conditions at control 
points on the surface. For vortex panels the equation to be satisfied is, 

ntv ntb 

2_ a ”j V + Y. D “j "p - ' v ■ v ■ w °. + < d v dx) i 

j-i J j-l iii 


and for body panels on the internal boundary it is, 


ntv 




ntb 

r 

j-i 


of. 



* o 


i 


The known perturbations from others singularities have been placed on the 
right hand side. Corresponding sets of equations may be written for 
symmetrical or antisymmetrical loading. 
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UNIT SOLUTION BOUNDARY CONDITIONS 


Several types of basic and unit boundary conditions are considered and 
can be classified as either symmetric or antisymmetric. Linearized theory 
allows the superposition of these basic unit solutions. The p, q and r 
rotary derivative boundary conditions are the result of placing the 
configuration at a * 0 J « 0 in a flow field rotating at one radian per 
second. 


Symmetric : 

la) Basic - vortex panels (dz /dx) - w 0 - w 0 - w 0 

C B r <t 


(dz^/dx) - surface slope due to twist and camber 

w 0 - normalwash induced by slender 
B body thickness and camber 

w 0 — normalwash induced by thickness 

r source panels 

w 0 - normalwash induced by body 

a source panels 


a - - ( U * n ) + U 
oo n 


e • n ) + U - - n + U 
x n x n 


lb) Basic - body panels (internal boundary) - <f> 0 - <f> 0 

T O 


4> o ** velocity potential induced by thickness 

r source panels 

4> 0 - velocity potential induced by body 

o source panels 
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2a) Unit alpha - vortex panels 


ttt cos 0 - w - w 

180 c ot n o. 

B a 


dihedral angle 


w 


w 


- normalwash induced by slender body 
*B at unit alpha 

* normalwash induced by body 
source panels at unit alpha 


a 


- ( U • n ) 

co 


180 ( V " > 


180 n z 


2b) Unit alpha - body panels 


- K 


velocity potential induced by 


3a) Unit q rotation - vortex panels 


a body at unit alpha 
2 


(x-x ) cos $ - w - w 


eg 


*B 


w = normalwash induced by slender body 
^B undergoing unit q rotation 


w 


= normalwash induced by body panels 
undergoing unit q rotation 


- ( u 


r> L 


( x “ x ^r.) ■ ( Z * Z ,J n 

eg z eg x 


3b) Unit q rotation - body panels 


velocity potential induced by body 
panels undergoing unit q rotation 


4a) Unit flap - vortex panel 


180 


4b) Unit flap - body panel 


1 . for flap panel 
0. for others 

0 a - 0 
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Antisymmetric : 


la) Unit beta - vortex panels 


180 Sln 6 c ' ‘ W 0 

B a 


0 “ dihedral angle 


w - normalwash induced by slender body 

P B at unit sideslip 

w - normalwash induced by body 

« source panels at unit sideslip 


a - - ( U • n ) 

ao 


180 ( ®y* n > 


7T 

180 n y 


lb) Unit beta - body panels 


- 4 > 




<t>„ = velocity potential induced by 

p a body 


2a) Unit p rotation - vortex panels 

2 


<y*y ) cos 8 + (z-z ) sin 8 - w - w 

L c & c eg c J p ] 


B Pa 


w = normalwash induced by slender body 
^B undergoing unit p rotation 

w = normalwash induced by body panels 
undergoing unit p rotation 


a - - ( U • n ) 

OO 7 


~T~ (y-y ) n - (z-z ) n 
b L eg z eg' y J 


2b) Unit p rotation - body panels - ^ 


velocity potential induced by body 
panels undergoing unit q rotation 


28 



3a) Unit r rotation - vortex panels 


2 

b 


(x-x ) sin 6 
eg c 



w 

r 


a 


a 


w = normalwash induced by slender body 
B undergoing unit r rotation 

w — normalwash induced by body panels 
a undergoing unit r rotation 

( U • n ) - - “““ \ (x-x ) n - (y-y ) n 1 
« b L eg y y eg x J 


3b) Unit r rotation - body panels 



a 



a 

“ velocity potential induced by body 
panels undergoing unit r rotation 


4a) Unit flap - vortex panel 


7T 

180 K 


k = 1. for flap panel 

k = 0 . for others 


4b) Unit flap - body panel 


0 a * 0 
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CONSTANT SOURCE AND CONSTANT VORTICITY PANEL INFLUENCE EQUATIONS 

The source finite elements have a discontinuity in normal velocity 
across the panel surface while the vortex finite elements have a 
discontinuity in the tangential velocity in a direction normal to the panel 
leading edge. The magnitude of the discontinuity, in each case, is constant 
over the panel area. In addition the vortex panels have a system of 
trailing vorticies extending undeflected to downstream infinity. 

A constant pressure or constant source panel with a quadrilateral shape 
can be constructed (figure 7) by adding or subtracting four semi - inf inite 

triangular shaped panels. These semi- inf inite triangles, each determined 
by a corner of the quadrilateral, can be assumed to induce a velocity 
perturbation everywhere in the flow. However, each corner represents only 
an integration limit, and all four corners must be included to make any 
sense. These perturbation potential expressions are derived in Appendix A. 

2>(x,y,z) = <Mx-x 1( y-y lf z,T 21 ) - tf(x-x 2 , y-y 2 , z,T 21 ) 

- *(x-x 3 , y-y 3 , z,T 43 ) + ^(x-x 4 , y-y 4 , z,T 43 ) 



Figure 7. Constant Pressure or Constant Source Panel Construction. 
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For one corner, having sides determined by y - 0 and x-Ty - 0: 

2 


2 2 2 2 2 2 2 
r = y + z , R - x + £ r , k 


r 1 p > 0 ' 

2 ' > 
2 p < 0 w 


2 2 2 2 2 
P - 1-M , B * T + p 


constant source panel 


^ s (x,y,z,T) 


ak 

4tt 


r i R+x 1 ^ BR+(Tx+/? y) 

I y 2 lo B + (x-Ty) - 2 lo S 5 — 


R-x 


B 


BR-(Tx+£ y) 


zR 


+ z tan 


ak 


l 


xy-Tr 
2 

BR+(Tx+£ y) 


} 


u s (x,y,z,T) = - — ~ 7 log 2 

4 n B BR-(Tx+/3 y) 


v*. ^ ! ri-x 1 j BR+(Tx+£ y) % 

v s (x,y,z,T) \ 7 log - T - 7 log — l 

A/ir V D_v II QD / r T'^,_i _0 ,t\ J 


a k 
4 tt 
ak 


w s (x ,y , z , T) 


tan 


4tt 


R+x 
R-x 
L zR 
xy-Tr 


BR- (Tx+£ y) 


constant vorticity panel 

kCp r 1 R+x 

-I Tz 7 log 

fi " r ^ R-x 


^ v (x,y,z,T) 


z B 


2 1 , BR+ ( Tx+/5 y) 


log 


B 


zR 


+ (x-Ty) tan 


BR- (Tx+yS y) 

2 i y 


T * ( 2_k ) [ Tz 7 l°g r + (x-Ty) tan - ] \ 
xy-Tr z J 


kCp j 
8 tt '• 


u v (x,y,z,T) 


zR 


tan 


xy-Tr 


( 2 -k) tan 


,i y 


} 


kCp , j zR zR 

v (x,y,z,T) “ - i T tan" o + 

v 


w v ( x , y , z , T ) 


87T 

kCp 

8tt 


2 T 2 
xy-Tr r 


i y zx > 
(2-k) [ T tan - - -j ] f 

7 r 


z r 


{ 


R+x 2 T i BR+ ( Tx+/3 y) yR 

log - B - 7 log j — + — 

R-x B BR- (Tx+yS y) r 


i 2 yx 

- (2-k) [ T 7 log r 


r J 


Only the real (not imaginary) , downstream, contributions are considered when 


M > 1 . 

GO 
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CONSTANT SOURCE AND CONSTANT DOUBLET PANEL INFLUENCE EQUATIONS 

Source and vortex panels used to represent body shapes may have an 
arbitrary quadrilateral shape, i.e. they need not have two streamwise edges. 
The influence equations may be written in the z - 0 plane, and a coordinate 
transformation used to obtain the perturbations of a panel having arbitrary 
orientation (see Appendix C) . A quadrilateral source panel of arbitrary 
shape can be constructed by combining quadrilaterals with streamwise 
parallel sides. 


(*i »Yi ,0) 



■* 

^ s ( x ‘ x l . 

y-yi . 

z 1 ) 

- ^ S ( X ' X 2 . 

y-y 2 . 

Z ,T 2 j 

+ 

^ S ( X - X 2. 

y-y 2 . 

z * ^3 2 ) 

- ^ S ( X - X 3. 

y-y 3 . 

2 » ^32 

- 

^ S ( X ' X 4. 

y-y 4 - 

z > ^3 4 ) 

+ ^ s (x-x 3f 

y-y 3 . 

2 »t 34 

+ 

^ S ( X - X 4 > 

y-y 4 . 


- * s (x-x lt 

y-yi. 

z , T 4 i 

= 

^ S ( X ‘ X 1. 

y-y x . 

z > i ) 

- 0 g (x-x lf 

y-yi. 

Z >^41 

+ 

^ S ( X - X 2> 

y-y 2 . 

Z » ^3 2 ) 

- ^ S ( X - X 2 > 

y-y 2 > 

Z , Ti 2 

+ 

^ S ( X ' X 3 - 

y-y 3 > 

2 >^ 43 ) 

- ^ S ( X - X 3. 

y-y 3 . 

z > T 2 3 

+ 

^ S ( X ' X 4- 

y-y 4 . 

2 , Tl 4 ) 

- ( X-X 4 > 

y-y 4 . 

2 ? T 3 4 


Therefore each corner consists of the difference between the perturbations 
induced by the two sweep angles. Therefore we can omit terms independent of 
T, since they will cancel when the two contributions are combined. 
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Therefore, ommiting terms independent of T, the perturbation velocities and 
perturbation potential for an arbitrary quadrilateral constant source panel 
are : 

Constant source panel 


ak r 1 x BR+(Tx+/? y) 

$ (x,y,z,T) - - — \ (x-Ty) - j log — 

4ir 1 B BR- (Tx+/3 y) 

2 

ak 1 l BR+(Tx+£ y) 

u (x,y,z,T) - - — - 7 log 2 — 

4 * B BR- (Tx+/J y) 


2 



ak 

T J BR+(Tx+/3 y) 

v s (x,y,z,T) - 

— 

“ 1 lo S 2 


4tt 

B BR - (Tx+y9 y) 


ak 

x zR 

w s (x,y , z , T) - 

- — 

tan 2 


47T 

xy-Tr 


A constant doublet panel is obtained by taking the z derivative of the 
constant source panel. 

^uk 1 zR 

^ d (x,y,z,T) - - — tan' y 

4tt xy-Tr 

2 

Mk 1 z (Tx+/3 y) 

u (x,y,z,T) -- — - 2 y - ! - 

4?r R [(x-Ty) + B z ] 

2 

/ik 1 Tz (Tx+jS y) 

V, (x,y,Z,T) “ - 2 2 2 

4?r R [ (x-Ty) + B z ] 


1 (x-Ty)(Tx+/y) 

w (x,y,z,T) - - — - j 2 ~T~ 

4?r R [ (x-Ty) + B z ] 
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Since the sweep angle could become infinite, we can write the above 
equations in a different form. First, the sweep angle can be written as 


T AX 

then 

A 2 

define B - 

2 

AY ’ 

AX 

where e.g. for T - T 12 

AX - 

x 2 - x x 

AY - 


Now for the source panel, using the previous definitions, we can write, 

{ 


^ s (x,y,z,T) 


a k 
4?r 


1 x BR+ (xAX+jflyAY) 
( x AY - y AX ) 7 2 log 2 

B BR - ( x AX + p y AY ) 


AY zR 


z tan 


xy AY - r /XX 
2 


} 


u s (x,y,z,T) 


v s (x ,y , z , T) 


w s ( x .y>z.T) 


ok AY j BR + ( x AX + /? y AY ) 

-a - J log a j 

4tt B BR - (xAX + /?yAY) 

A 2 

ok AX x BR + ( X AX + p y AY ) 

* 2 1 °S 2 


47r B 


ak 

4tt 


— tan 


BR - ( x AX + y AY ) 
t AY zR 

2 

xy AY - r AX 


The constant doublet panel is now. 


^ d (x,y,z,T) 


u d ( x > y > z , t ) 


v d (x,y,z,T) 


w d (x,y,z,T) 


Mk AY zR 

— tan 2 — 

4tt xy AY - r AX 


2 

pk 1 z AY ( x AX + $ y AY ) 


4tt 

R 

[ ( 

2 

x AY - y AX ) 

a 2 2 2 

+ B z AY ] 

/ik 

1 


z AX ( x AX + 

/y ay ) 

47T 

R 

[ ( 

2 

x AY - y AX ) 

A 2 2 2 
+ B z AY ] 

^xk 

1 

( x AY - y AX ) ( x 

2 

AX + $ y AY ) 

47T 

R 

[ ( 

x AY - y AX ) 2 

a 2 2 2 

+ B z AY ] 
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LINEARLY VARYING SOURCE PANEL INFLUENCE EQUATIONS 


In supersonic flow constant source panels having a sonic edge have a 
real singularity along an extension of this edge. The singularity occurs 
because : 


lim 

(x-x x ) - T(y-yj) -+ 0 


2 2 2 
e “ (T +/})-* 0 



i cR 1 +[T(x-x 1 )+£ (y-yj) ] 

7 lo 6 i 

£ Ri - [T( x-x 1 )+^ (y-y x ) ] 


i €R 2 +[t(x-x 2 )+0 (y-y 2 ) ] 

7 lo s 5 

€R 2 - [T(x-x 2 )+/3 (y-y 2 ) ] 


} 


co 



Control points which are near the extension of this edge will have large u 
and v velocities induced upon them. The singularity can be eliminated by 
using panels which have a source distribution which varies linearly in the 
chordwise direction. The resulting continuous source distribution 
eliminates the singularities. The linearly varying source panel influence 
equations can be found by integrating the constant source panel influence 
with respect to x. 


u 


1 o 


V 


1 0 


-k 

r j R+X 


i i 

< 

y 7 i°s 

+ 

(x-Ty) - 2 log 


R-x 


B 

_k i 

r 1 

R+x 

1 

- 

(x-Ty) 2 log 


- T(x-Ty) - 

2n ' 


R-x 

B 


2 

BR+(Tx+£ y) 


2 ' ^ 
BR-(Tx+0 y) 

2 

! BR+(Tx+/3 y) 

7 log 2 — 

BR-(Tx+/9 y) 


tan 


x zR 
xy-Tr 


} 



' k r 

1 

R+x 

W 10 

= — { tz 

7 lo S 





R-x 


- k r 

1 

R+x 

^10 

= - { 12 

7 1°S 



2n ^ 


R-x 


_i zR a 

Tz tan 2 ' R r 

xy-Tr ' 

2 

2 1 i BR+ ( Tx+/3 y) 

- B Z - 2 log 2 — 

B BR-(Tx+0 y) 

2 

2 1 i BR+(Tx+/9 y) 

- B z — 7 log 2 — 

B BR-(Tx+/9 y) 


(x-Ty) 

(x-Ty) 


! zR 

tan 

xy-Tr 

t zR 

tan 

xy-Tr 


} 

} 
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These velocity components satisfy the same criteria as the velocity 
components for the constant source panels except that the source strength is 
proportional to (x-Ty) . The source panel finite elements are constructed 
with the following properties. 

1. All panel leading and trailing edges are at constant (x/c) , side 
edges are at constant y. 

2. Each source finite element is composed of a pair of chordwise 
adjacent panels. 


3. The source strength varies linearly with chord measured from the 
leading edge of a panel pair, i.e. the maximum value of the source 
strength is proportional to the local chord and attains this 
maximum on the panel edge joining the panel pair. 


y 

x 


i 


5 


\ 

\ 

\ 

\ 

\ \ 


\ 

\ 


A = (x/c) - 

(x/c) - 

(X/c) - 

(x/c) 

\ 

\ 


3 1 3 

1 

4 


\ 

\ 






\ 

2 

A - (x/c) - 

(x/c) - 

(X/c) - 

(x/c) 

\ 


S3 5 

3 

6 



\ 


2 


4 


\ 


a 


4 

6 


The perturbation velocities induced by this panel pair are composed of 
contributions from six corners. 


u(x,y ,z) - <j/A 31 [u 10 (x-x 1 , 

y-yi . 

z > 

Ti) 

- u 10 (x-x 3 , 

y-y 3 . 

z, 

t 3 ) 

-u 10 (x-x 2 , 
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If there are N panels in the chordwise direction there will be N-l 
singularities or unknown source strengths associated with them. The linear 
variation in the source distribution means the value of dz/dx must be zero 
at the leading and trailing edges of each span station. This may be an 
undesired restriction and therefore the use of linearly varying source 
panels is optional. 
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EDGE EFFECTS 


The low pressure created by high velocities around a surface subsonic 
leading edge results in a suction force. As the edge becomes thinner or the 
angle of attack increases, the flow deviates from potential conditions 
resulting in a progressive loss of theoretical suction and an increase in 

5 

drag. Generalizing a concept due to Polhamus , it is assumed that the 
leading edge vortex created by the detached flow in effect rotates the lost 
suction force perpendicular to the local surface. 

In order to implement this philosophy, a method of determining the 
spanwise variation of potential suction was developed using linear thin wing 

theory and involves finding the coefficient of the l/7x term in the 
chordwise net pressure distribution. The analysis is applicable to multiple 
surface problems of arbitrary planform in the presence of bodies at any Mach 
number. If the chordwise net pressure distribution on a thin wing at any 
given span station is expanded in a series 


N 

AC p = A„ cotHp) + A n sin(n^) 

n-1 

£ = x/c - —■ (1-cos <j > ) - sin 


( 12 ) 


it is shown in appendix B that the leading edge nondimensionalized suction 
force per unit length is 


_ , . ATHRUST 

C s (y) - — iy-,: 


_7t_ 

8 


2 2 
J T + fi 


where T - tan A T „ ; B - 1-M 

L . h « 


and c is the local chord. 


(13) 


Only the first term in equation 12 contributes to the thrust, since it 
is the only one which is infinite at the leading edge. If the chordwise 

pressures are known at M points along the chord, the coefficients A are 

t n 

obtained by fitting a least square error curve described by N terms of the 
series, through the points, where N < M. The pressure distribution is 
obtained using constant pressure panel analysis. 
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The method used to compute the suction force at surface tips is similar 
to that for the leading edge. By using the irrotational property of the 
flow, it is shown in appendix B that the tip suction force is 


C s ( ° " AF „ /(q« c T Ax) ■ (7r/32)tc lvG /(c T ,? MAX )lc n 0 f/f(*)dx] 2 (14) 

where 
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tip surface lateral surface dimension 
faction of chord 

c n (,)[c(,)/c AVG I, MSX [ 0,;^ - , 2 )/2f 1/2 
is local section normal force coefficient 


and as ij -* the net pressure coefficient is assumed to be of the form 


f 2 ) 1/2 

AC p (f “ { [ 1 ' (,?/, 7 MAX ) 1/2 } (C AVG /C(,?)) C n 0 f(O 


The sectional leading edge suction attained in the real flow, C (y) , is 


estimated by 


C s (y) - K s (y) C s (y) 
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and 


M - - J2 C" 1 [(1+C 2 ) 1/2 - 1] 1/2 
e p p J 

5 p - C p,LIM ; ~ (1 ' M n > V2 

C p,LIM"' 2/(7M n )[(R n X lo6 >/(R* 10 ' 6 + 10 (4 ' 3M n ) ) j ( - 05 + - 3 5 ( 1 -M r ) ) 

R n “ R - (c n /a> cos(A LE ) 
c 

2 

C — C (c/c ) / cos (A__) 
s,n s v ' n 7/ v LE' 

M = M cos (A t _) 
n oo v LE 

The chord of the normal section, c , is defined so as to place the 
maximum thickness, t , at the mid-chord as indicated in figure 8. The 
associated leading edge radius is designated by r 


Leading 



Potential tip suction is assumed to be fully rotated as a result of 
vortex formation in the present analysis. 
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JET FLAP 


A completely linearized approach is used based on the assumptions of 
thin airfoil theory. The flow is assumed to be inviscid and irrotational 
and all entrainment effects are neglected. The jet is represented by an 
infinitesimally thin sheet having zero mass flow but finite momentum per 
unit of span. This sheet is assumed to extend from the trailing edge of the 
surface back to infinity. (In practice one or two chord lengths is 
sufficient) . The effects of transverse momentum and the deflection of the 
jet sheet are neglected. 

Since both the planform and jet can maintain a pressure discontinuity, 
they are both represented by a system of quadrilateral panels having 
continuous distributions of vorticity. The strengths of these constant 
pressure vortex panels are determined by solving a set of linear 
simultaneous equations which satisfy the downwash boundary conditions at a 
set of control points on the planform and jet. 

The boundary condition on the planform is the previously described flow 
tangency condition. The pressure difference across the jet causes a change 
in the direction of the jet momentum. The equation relating these 
quantities forms the boundary condition on the jet and can be derived by 
considering a jet segment of unit depth. 



The mass rate of flow through the jet is m and the velocity is V. If we 
assume a pressure difference of AP across the jet, then from the momentum 
theorem applied to the differential element, we write 


or 


mVjA0 = RAP t\<f> 


AP = mV. 

J 


/R 


where R is the radius of curvature of the jet. 


40 



The reaction of the jet on the flow external to the jet is 

F - RAP A <f> 

A vortex of strength 7 per unit length along the jet would produce a 
reaction of 

F “ p U 7 R A<f> 

00 


Hence, equating these two forces, we calculate the action of the jet on the 
flow external to the jet by replacing the jet with a running vortex strength 
7 given by 

7 - ®V /(/>UR) 

For a nearly horizontal jet with a large radius of curvature 

2 2 

1/R ~ dz/dx = dw/dx 

where w is local downwash velocity (nondimensionalized with respect to U ) . 


Then 


or 





C 




pU dx 

00 


c M (y) c (y) w(x,y) - AC p (x,y) - 0 


which is the boundary condition written for a three dimensional jet flap. 

To apply the jet flap boundary condition to control point i, the above 
equation is integrated between adjacent control points in the streamwise 
direction. 


c/y) C(y) 


x 

c . 
1 



w(x,y) 


dx 



I 


AC p (x,y) 


X 

c 


i -1 


dx 


= 0 


(15) 


The control point is located at 87.5 percent of each panel chord. To 
simplify the second integral in equation (15), the assumption is made that 
the control point is exactly at the panel trailing edge. The effects of 
such an assumption have been shown to be negligible. Equation (15) 
evaluated from the leading to the trailing edge of panel i yields the 
following relation: 


C C 




Vi 1 


ac: Ax. ** 0 

p. 1 
1 


(16) 
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The downwash at each control point is written in terms of the N net 
pressures on the quadrilateral panels: 


w. 
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Equation (16) is then written 
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where 5., is the Kroneker delta, 
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P . 
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For a flap panel adjacent to the jet exit, equation (16) must include 
any jet deflection angle relative to the surface trailing edge. 


or 
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is the jet deflection angle. 
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The complete set of linear simultaneous equations for both the surface 
and the jet flap is then written 



i - 1,N 


where 



A. . for i on the 


C C [A, . -A, 14 ] - 
/i ij i-lj J 


surface 


8 . . Ax , 

i 


for i on the jet 


(17) 


for i on the surface 

C^C<5^ for i on the jet adjacent to the exit ► 

0 for i elsewhere on the jet 

Both symmetric and antisymmetric jet deflections are considered. 
Thus, after calculating the influence matrices and boundary conditions in the 
usual manner, the appropriate rows are modified and combined to produce a 
linear symmetric or antisymmetric system as described by equation (17). 
Because of the rotational quality of the flow fields, the p, q and r rotary 
derivative calculations are generally not valid for jet flap configurations. 




INLETS 

A jet boundary between two flows with different total energies is 
characterized by a discontinuity in tangential velocity, but a continuous 
value of Cp . This flow can be replaced by a flow with the same total 
energy everywhere, but now having a discontinuity in Cp across the jet 
boundary, instead of a discontinuity in total energy. The jet boundary will 
be represented by a vortex sheet having the same discontinuity in tangential 
velocity, and the velocities will be the same everywhere in the two flows. 
If the jet is such that the perturbation velocities are small compared with 
the free stream velocity, i.e. u,v,w « U^, this jet boundary can be 

simulated by constant ACp chord plane panels. It will be shown that, to 
first order, this value of ACp is constant on the entire jet boundary. Let 
u be the x component of the perturbation velocity. Then, assuming energy 
addition to the jet flow, write the energy equation across the jet boundary: 
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To first order, 
therefore 
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to first order the energy equation becomes, 
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00 J 
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+ 00 00 
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1 
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2 12 

(2-M ) TP U 
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const 


where ACp is in the direction of the normal pointing into the jet flow 

The inlet is simulated by specifying an average mass flux over a set of 
field points within the inlet region. The nondimens ional mass flux per unit 
area in the x direction can be written: 


£ U 
p U 

oo oo 


(1 - M u) (1 + u) 


[ 1 + (1-MJu ] 


The value of this expression can be calculated, on each field point, for a 
unit value of ACp across each panel of the configuration and jet boundary. 
Therefore a linear equation can be written to constrain the value of the 
average of the mass flux, and therefore the inlet mass flux, on a given set 
of field points. The additional unknown required for the additional 
equation is supplied by the (constant) value of ACp across the jet boundary. 
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AERODYNAMIC CHARACTERISTICS 


Longitudinal and lateral - directional forces and moments due to 
thickness , ^twist and camber, pitch, sideslip, and the dimensionless rotary 

velocities p, q, and r are obtained from surface pressure integrations of 
the various configuration components. 

Slender Bodies 


The pressure coefficient, to an approximation consistent with slender 
body theory, is 
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b/2 J 




(18) 


Paneled Bodies 

For paneled bodies a surface differentiation of the perturbation 
potential is used to obtain the perturbation velocity components tangential 
to the surface. The velocity normal to the surface is obtained from the 
imposed boundary condition. A formula for the pressure coefficient can be 
derived using the energy equation. Although the perturbation velocities 
were obtained using a linear equation, on body surfaces, which may be quite 
thick, a better approximation to the pressure coefficient is obtained using 
a nonlinear formula. Assuming the freestream is at an angle of attack a and 
an angle of sideslip /?, and the perturbation velocities are nondimensional , 
we can write for the freestream, 
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For the energy equation we can write, 
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'Y - 1 2 

1 - M 

2 00 


2 ( u cos a cos f} - v cos a sin p + w sin a ) 
w 2 ] 


2 2 2 
+ U + V + 


Since the existence of a velocity potential assumes isentropic flow, we can 
use the isentropic relationship between pressure and temperature. 
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for small values of a, and this becomes 
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Using the above expression the isentropic pressure formula for c is • 
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For small values of the quantity 8, 
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we can expand the exponentiation in an infinite series in 5. 
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Retaining only terms up to order 2 in a, p, u, v, w, we have, 


2 f 2222 

2 (u-/Jv + aw) + (1-M^) u + v + w 


or 
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Where, to first order, the freestream is represented by, 
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Planar Components 


Surface pressure distributions are calculated for planar components 
using the first-order linearized form 


-2u/U - -2 [(u/U) ind ± C 


/ 4 ] 


NET 


The ± signs refer to the upper and lower surfaces respectively. The term 
consists of the velocities induced by the isolated bodies and other 


vortex and source panels. These velocities are obtained by multiplying the 
(u/U) influence matrices by the appropriate panel strengths. The C /4 

P NET 

term accounts for the u/U perturbation velocity induced by the local 
distribution of vorticity and changes sign from upper to lower surface. The 
total (u/U) and C values are the result of taking linear combinations of 
NET 

all the basic and unit solutions. 


FORCES AND MOMENTS 
Slender Bodies 


The forces and moments are obtained from the surface integrations. Let ds 
be an element of arc length at a given x section, and let x be 
nondimensionalized with respect to the body length L. First performing line 
integrals at each section and then integrating over x gives, 
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In terras of these expressions, the commonly used aerodynamic 
coefficients are 
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where L is the body length and c, b and S are configuration reference 

KfcjT 

chord, span and area, respectively. 

Cross - coupling between the pitch, sideslip, and rotary motions through 
the product and quadratic terms in equation (18) is neglected. 


Paneled Bodies 


A surface differentiation done on the perturbation potential for each 
of the basic and unit solutions to obtain the velocity components tangential 
to the surface. The velocity normal to the surface is obtained from the 
imposed boundary condition. These unit solution surface velocities are 
combined to obtain the resultant pressure coefficient. To obtain the 
section forces and moments, component forces and moments, and configuration 
forces and moments, a surface integration is done. Each computed pressure 
coefficient is multiplied by the panel area and the proper component of the 
surface normal, and the result is summed over all of the body panels. 

Planar Components 

The perturbation velocities for each of the basic and unit solutions 
are combined to give the resultant pressure coefficient. The net pressures 
and pressure coefficients are then integrated numerically to give the 
section forces and moments, component forces and moments, and configuration 
forces and moments. 

Since the vortex panels have a constant pressure distribution, a block 
integration scheme is employed. With the exception of drag, these basic and 
unit force and moment coefficients are combined in a linear manner to 
produce the aerodynamic characteristics for any desired flight condition. 
Since drag varies in a parabolic manner, it must be considered on a point by 
point basis as defined in a later section. 

The longitudinal normal force distribution on the bodies is calculated for 
each solution. The load distribution on the interference shell portion of 
the body is given by integrating over all vortex panels at a given 
longitudinal station. 

Normal Force 


where N is the number of panels around the shell, L is the length of the 
body, Ax is the length of the interference shell segment, is the panel 

area, and C x * 2 for a centerline body or C t - 1 for an off centerline body. 

This carryover load distribution is added to the previously calculated 
isolated body longitudinal load distribution. 
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The section characteristics of planar components are determined by a 
chordwise summation of panel data at each span station and are given by the 
following equations: 


Local Normal Force: 
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Weighted Normal Force: 
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Weighted Lift Force: 
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where Nc is the number of chordwise panels, 6 is the section dihedral angle 

and As is the width of the span station and is given by 

2 2 1/2 
As - [Ay + Az ] ' 


The section characteristics due to the reaction of a jet flap are 
calculated by taking the appropriate component of the reaction force. 


Reaction Normal Force: 
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where £. is the total deflection angle of the jet. 


Reaction Lift Force: 
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Component forces and moments including edge vortex effects are given by 
the following equations : 
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Lift: 


C L S REF - F » 


N NS 

} \ ET A l“ S < # l> + } (1 - K s, )C s, C j iS j n i 
1-1 " tl i j -1 J J 


C <k A(x/C T ) k 


Side Force: 


Vref - F * 2 C P NET v in <V • 2 (1 V % 

- Ktii. . i J J J 


i-l i j- 

nc t 

• C T ) C s v 4(X / C I>k T Y„ 
k-1 k k 


Rolling Moment: 


C i bS REF “' F 2 ) C P A i [(y.-y CG )cos(fl i )+ (z.- 2 CG )sin((? i )] 


NET. 
i-l l 


NS 

■ } <1 - K s. ) c s. C j As j [( yLE-yCG> n L. + (z LE- 2 CG )n Y. 1 

j_l J J J J J 


C T 2 C S. (Ax/C T ) k ^ y T* y CG^ T L. + ^ Z t' Z CG )T Y J 
t-i k ^ k 


Pitching Moment: 


V S REF - F » 2 C P M , T A i (X i- X CG )COS(<? i ) 

# a. 1 NL 1 « 

i=i i 


(l-K ) c s e AsJ( x le - x cg) Q l.' C T / C s. A(x/C T ) k T Y, 
J J J J J k -l k k 


Yawing Moment: 
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where N is the number of vortex panels on half of a symmetrical component 
(or total for an asymmetrical component) and F x , F 2 are given by 
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symmetric loading: 


antisymmetric loading: 
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f 2 

Fi 
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* 1 asymmetric geometry 

— 2 symmetric geometry 

— 1 asymmetric geometry 

— 0 symmetric geometry 

= 1 asymmetric geometry 

* 0 symmetric geometry 

** 1 asymmetric geometry 

— 2 symmetric geometry 


For the leading and side edge vortex terms, Ns is the total number of 

spanwise panels for both component halves, NC T is the number of tip 

chordwise panels, is the axial location of tip vortex center of 

c ♦ P • 


pressure , As ” As./l + T and the rotation factors Q and T are derived in 
appendix B and defined below. 


Leading Edge Vortex Rotation: 

” “ s * n a ( cos A cos 8) + cos a (cos 9 sin 8 - sin 0 sin A cos 8) 

+ A 0 / | A 0 | [ sin a (cos A sin 6) + cos a (cos 9 cos 8 + sin 9 sin A sin 8)] 

Qy =* cos 9 sin A cos 8 + sin 9 sin 8 

+ A 0 /|A 0 |[sin 9 cos 8 - cos 9 sin A sin 8] 

where 8 is the slope angle of the leading edge camber line and the sign of 
coefficient A 0 (from equation 13) is used to determine the direction of 

vortex rotation. 

Side Edge Vortex Rotation: 

Ty ~ — cos a sin 9 + C / | C | (sin a sin 8 + cos a cos 9 cos 8) 

n o “o 

Ty - ± cos 9 + /|C^ I sin 9 cos 6 


where S is the slope angle of the tip camber line, ± is plus for the left 

side and negative for the right side of the configuration and the sign of 

coefficient C (from equation 14) is used to determine the direction of 
"o 

vortex rotation. 
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The x-coordinate of the center of pressure is given by 
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For interference shell components, the total forces and moments of the 
corresponding isolated body are added to those of the shell. 


Jet reaction forces and moments are obtained from a spanwise summation 
of the jet flap section characteristics: 
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where N is the number of spanwise jet flap stations and F 2 and F 2 are as 
previously defined. 

The forces and moments for the complete configuration are obtained by 
summing those of the individual components . 
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DRAG ANALYSIS 


Estimation of configuration aerodynamic efficiency requires the 
calculation of drag. The analysis separates the computation into skin 
friction and pressure drag components that are assumed to be independent of 
each other. The following form is considered and produces nonparabolic 
polars as a result of the incorporation of edge force considerations . 


°D- C D. +C D + C n +C n 

viscous wave base lift 


The specific techniques used for the various drag evaluations are discussed 
below. 


SKIN FRICTION 


Several well established semiempirical techniques for the evaluation of 
adiabatic laminar and turbulent flat plate skin friction at incompressible 
and compressible speeds are used to estimate the viscous drag of advanced 
aircraft using a component buildup approach. A specified transition point 
calculation option is provided in conjunction with a matching of the 
momentum thickness to link the two boundary layer states. For the 
turbulent condition, the increase in drag due to distributed surface 
roughness is treated using uniformly distributed sand grain results. 
Component thickness effects are approximated using experimental data 
correlations for two-dimensional airfoil sections and bodies of revolution. 

Considerations such as separation, component interference, and discrete 
protuberances (e.g. antennas, drains, aft facing steps, etc.) must be 
accounted for separately. 

In the following, a discussion is presented for a single component 
evaluation in order to simplify writing of the equations and eliminate 
multiple subscripting. The total result is obtained by a surface area 
weighted summation of the various component analyses as described on page 
44 . 


Laminar /Trans it ion 


A specified transition option is provided in the program. The 
principal function of the calculation is to provide the conditions required 
to initialize the turbulent solution. In particular, the transition point 
length and momentum thickness Reynolds numbers are required. 


R 

X TRAN 


R(x, 


TRAN 


/L)L 


R, - 0.664 7r C* 

TRAN X TRAN 
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where 


c* = (Z/^xiyx*) 

T/T^ - 1+0.72 [(T /T^) - 1] 

T r /T a - 1 + y?r (7-D/2 M 2 - 1 + 0.851( 7 -l)/2 M 2 
A* - 2.270 x 10' 8 T 3//2 /(T+198.6) lb-sec/ft 2 


This solution is based on the laminar Blasius result (8, chapter VII) in 

q 

conjunction with Eckert's compressibility transformation . This option 
permits an assessment of the reduction in skin friction drag if laminar flow 
can be maintained for the specified extent. It does not establish the 
liklihood that such a condition will be realized in practice or to what 
extent . 


Turbulent 


Smooth and distributed rough surface options have been provided in the 
analysis. In either case, the solution is initialized by matching the 
momentum thickness at the transition point produced by the laminar solution. 
That is, an effective origin (commonly referred to as a virtual origin) is 
established for the turbulent analysis. 

For the hydraulically smooth case 


C- R A 
F Ax 


2 R, 


TRAN 


R 


Ax 


C F R Ax / C F 


Ax 

i 

R. 


R. / R 

Ax 


L ' ^TRAN + 


(R)(J) 


solve for C' using, 
r 


C from equation (19) 

r 

for known R, 

F Ax 


0.242[sin' 1 a + sin'W [ ( 7 -l)/2M 2 C^] 1/2 = log^C^) - co log 10 (T^) 
then, 

C F " 29 x^o /L " <2« TE /i)(i/L) - Cp(i/L) 


( 19 ) 
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where 


TRAN 


R _1 R 


TRAN 


(2A 2 - B)[B 2 + 4A 2 ]' 1//2 


B[B 2 + 4A 2 ]" 1 / 2 


(7-D/2 (TyTp 


[1 + (7-D/2 M co ](T oo /T r )- 1 


0.88 


U) 


0.76 


The compressible turbulent flat plate method used here is that proposed 


1 o 


by Van Driest based on the Von Karman mixing length hypothesis in 
conjunction with the Squire-Young formulation for profile drag (8, chapter 
XXIV) as applied to a flat plate. 


For the distributed rough case 


Ax 


TRAN 


[1.89 + 1.62 log 10 (Ax/K s )]' 2 ‘ 5 [l + r( 7 -l)/2 M 2 ]' 1 


Ax 


2 C F 1 R 0 

r rRAN 


Ax 


R" 1 R 


Ax . - 
l+l 


L - X TRAN + AX 


[1.89 + 1.62 lo gl 0 (i/K^) ] 2 ‘ 5 [1 


r(7-l)/2 M 2 ]' 1 


C^. (i/ L) 


MAX [C , C ] 

SMOOTH ROUGH 
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The turbulent flat plate method used here is that of Schlicting (8, 
chapter XXI) which is based on a transposition of Nikuradse's densely packed 
sand grain roughened pipe data. The effect of compressibility is due to the 

reduction in density at the wall as proposed by Goddard/ 1 The selection of 
the equivalent sand grain roughness for a given manufacturing surface finish 

is made with the aid of Table II which was taken from Clutter. 


TABLE II 


Type of Surface 


Equivalent Sand Roughness 
(inches) 


Aerodynamically smooth 
Polished metal 
Natural sheet metal 

Smooth matte paint, carefully applied 
Standard camouflage paint, average application 
Camouflage paint, mass -production spray 
Dip -galvanized metal surface 
Natural surface of cast iron 


0 

0.02 - 0.08 x 10 
0.16 x 10' 3 
0.25 X 10' 3 
0.40 x 10' 3 
1.20 x 10" 3 
6 x 10 3 
10 x 10' 3 


Thickness Corrections 


The foregoing evaluations produce an estimate of the shearing forces on 
a flat plate (at zero angle of attack) for a variety of conditions. As an 
actual aircraft has a finite thickness, an estimate of pressure gradient 
effects on skin friction and boundary layer displacement pressure drag 
losses is required. A common procedure for accomplishing this and the one 
which will be used here is based on non-lifting experimental correlations 
for symmetric two-dimensional airfoils and axisymmetric bodies. The 
following relations derived by Horner (13, chapter VI) are used, 
respectively. 


K - C d /(2C p ) - 1 + K, (t/c) + 60 (t/c) 4 

- V C D "1 + 1.5 (d/L) 3/2 + 7 (d/L) 3 

F 

Horner recommends K : - 2 for airfoils with maximum thickness at 30% 

chord and K x - 1.2 for NACA 64 and 65 series airfoils. In this regard, the 

best information available to an analyst for his particular contour should 
be used. This is especially true for modern high performance shapes such as 
the supercritical airfoil. 
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Total Viscous Drag 


The aircraft total viscous-drag coefficient is estimated by a sum of 
the preceding analysis over all components (i.e. wing, fuselage, vertical 
tail, etc.)* That is 


C 


D . 

viscous 


1 c F.<yw t J 

j-i J 


The component length used in the calculation of the skin friction 
coefficient is the local chord for planar component segments and the 
physical length for bodies and nacelles. 


BASE DRAG 

Blunt base increments are estimated at subsonic and supersonic speeds 
by 

Cp BASE S ' ASE ' SrEF 
0.139 + 0.419 (M - 0.161) 2 ; M < 1 

CO 03 

M' 2 - 0.57 M‘ 4 ; M > 1 

CO 00 * co~ 

The expressions for the base pressure coefficient are derived from 
correlation of flight test results for the X-15, various lifting bodies, and 
the space shuttle. Power effects are treated as reductions in base area in 
the present analysis. 


where 


AC 


BASE 


-C T 


BASE 


'BASE 
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POTENTIAL DRAG 


One hundred percent suction drag due to lift and supersonic wave drag 
due to thickness can be evaluated by integration of the momentum flux 
through a large circular cylinder centered on the x-axis and whose radius 
approaches infinity (figure 9) . 


z 



Figure 9. Integration of Momentum Flux Through Large Circular Cylinder. 


The resulting expression for the total pressure drag is 


C D S REF 


-2 ff $ x $ r dA 2 + //($* + $ x )dA 


C D S REF 
w 


C D V S REF 


The first term represents the wave drag due to momentum losses thru the 
side of the cylinder caused by standing pressure waves at supersonic speeds. 
The second term represents the vortex drag which arises from the kinetic 
energy left behind in the Trefftz plane by the system of trailing vorticies. 


Vortex Drag 

The vortex drag may be computed when the distribution of trailing 
vorticity in the Trefftz plane is known. The assumptions of linearized thin 
wing theory result in a vortex sheet which extends directly downstream of 
all lifting surfaces. By changing a surface integral for kinetic energy to 
a line integral over the vortex sheet in the Trefftz plane the following 
expressions for lift and drag result. 
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where 


C L “ (C AVG /S REF ) l C n (n)c° s HD dr, 

C D - <W S REF>J[ 

V (j 

C vortex sheet contour 

C weighted section normal force coefficient C (c/c AT7/ J 
n & n AVG 

w^ asymptotic normal velocity on the vortex sheet 
r) vortex sheet branch coordinate 

$ inclination of vortex sheet with respect to y-axis 
The analysis computes the normal velocity on the vortex sheet, w^ 

i 

by assuming the vortex sheet is composed of finite trailing horseshoe 
vorticies whose strength is proportional to the local section c^(s) . The 

normal velocity is computed at a control point located midway between the 
trailing vortex segments (figure 10) . 



Figure 10. Trefftz Plane Vortex Wake Nomenclature. 
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V ■ C AVG r j/ (2 * r > e * 


n % (y i+l ,Z i+l ) 


\ t / 


\ 


\ 


r - [c - c ]/2 

n. n. ' 

J J 


<y j :i j' 5 


/ 


/ 

/ 


, (y c.- Z c.> 

/ 11 


/ 


/ 


/ 


< y i’ z i> 


r - [<y c .- yj)e y + Zj)e z ] 


‘j, ■ t-< z ^ • z () e „ + <y„ -y,)e 1/r 


iS i - <y i + r y i )e y + 


Ay . e + Az . e 
i y i z 


c i J y c i J ZJ 


-r - - -♦> -f 

n =*= As. [-Az.e + Ay.e 
x l y 2 


w = V*n - c T 
<*>. AVG 

i 


Therefore 


where 


/(27rr A Si )|z c - Zj ^ + <y c - y^AyJ I\ 

C - ) w C A> 7 . 

D / oo. n. l 

.11 
l 

w - ) a, .r. 

®- L ij J 


a, .r. 
ij J 


Wave Drag 

The integral for wave drag 


C S = - 2 f f <f> <f> rdxdtf 
D, t REF “ j x r 

W A 2 


may be simplified by allowing the cylindrical surface of integration to 
recede infinitely far from the disturbance. Under these conditions, the 
spatial singularity simulations can be reduced to a series of one- 
dimensional distributions. The basis for this reduction is the finding by 
Hayes (14) that the potential and the gradients of interest induced by a 
singularity along an arbitrary trace on a distant control surface, say PP' 
of figure 11 (or alternately described by the cylindrical angle 8), is 
invariant to a finite translation along the surface of a hyperboloid 
emanating from the trace and passing through the singularity. As the apex 
the hyperboloid is a great distance away, the aforementioned movement is 
along a surface which is essentially plane; it will be henceforth referred 
to as an "oblique plane". Since a singularity is a solution of a linear 
differential equation, all singular solutions which lie on the surface of 
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Figure 11. Distant Control Surface Geometry 



the same hyperboloid (oblique plane) may thus be grouped to form a single 
equivalent point singularity whose strength is equal to the algebraic sum of 
the individual strengths and which induces the same potential (momentum) 
along the trace as the group of individual singularities. 

This finding provides the basic technique for reducing a general 
spatial distribution of singularities to a series of equivalent lineal 
distributions. This is accomplished by surveying the three-dimensional 
distribution longitudinally at a series of fixed cylindrical angles, 9 . At 
each angle, the survey produces an equivalent lineal distribution by 
systematically cutting the spatial distribution at a series of longitudinal 
stations along its length. At each cut, the group of intercepted 
singularities is collapsed along the "oblique plane" to form one of the 
equivalent point singularities comprising the lineal distribution. 

The far- field expression for the wave drag of a general system of lift, 
and side force elements is 

oo co 

^REF ** (4*U) / / S *^ ( e x , 0 )h^( e 2 , 0 ) ln| e x - e 2 | deide 2 d 0 

where 

h (e, 0 ) - f(e, 0 ) - g (e, 0 ) sin 9 - g (e, 0 ) cos 9 

“ z y 

is the equivalent lineal singularity strength at 

the cylindrical angle 9 

f (e, 0 ) =* equivalent source strength per unit length 
P 1 U “ equivalent lifting element strength per unit length 

P U gy(e, 0 ) = equivalent side force strength per unit length 

These strengths are deduced from the three dimensional singularity 
distributions by application of the superposition principle along 
equipotential surfaces. For a distant observer such surfaces are planar in 
the vicinity of the singularity configuration. The individual singularity 
strengths are related to the object under consideration by the requirement 
of flow tangency at the solid boundary. Lomax (15) derived the following 
approximate expressions between the equivalent singularity strengths and a 
slender lifting object. 

f (€,*) - u a/de[k (€,*)] 

g z <«.*> = <0U)/2 J c C p dy 

gy(e.tf) = (£U)/2 J c Cpdz 

where (see figure 12 ) 

A(e,0) - is the Y-Z projection of the obliquely cut 

cross - sectional area 

c is the contour around the surface in the oblique cut 
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Utilizing the singularity strength expressions derived by Lomax, the 
following expression for wave resistance based on the far field theory of 
Hayes is obtained 


C D S REF 
w 


2tt L(0) L (0), 2 2 

— — i — 2 J f f \d /de^kie^d)] 

(4*) 0 -L(8) -L(0) 1 


- /3/2 d/de^sind / C (e 1 ,6) dy + cos 0 J C 

c P c P * 

/de 2 [A(c 2 » 9) ] - /9/2 3/d£ 2 [sin0 J C (e 2 ,#)dy + cos0 J C (e 2 ,0)dz]j 


In | « j - e 2 | de j de 2 d0 


In order to facilitate subsequent discussion, the above result is 
manipulated into the following form 


2 7T 1 1 

C D / J J d / de i [ A „ ( € i . ]3 /ae2[A o (e 2 ^)]ln|£ 1 -£ 2 |de 1 d€ 2 d<? 

W ooo e e 


where the effective A is given by, 


( 20 ) 


A(e,0) “ A(e,0)- p/2 J f C (e,0)[sin ^ dy + cos ^ dz] de 


A requirement for this transformation is that 


A e (CM) = A e (L > 0) " 0 

In accordance with equation 20, the wave drag of a configuration is the 
average of the wave drag of a series of equivalent bodies of revolution. 
The drag of each of these bodies is calculated from a knowledge of its 
longitudinal distribution of normal cross-sectional area. For each 
equivalent body, these areas are defined to be the frontal projection of the 
areas and the accumulation of pressure force in the theta direction 
intercepted on the original configuration by a system of parallel oblique 
planes each inclined at the given Mach angle. The common polar angle (9) of 
the system identifies the equivalent body under consideration. 

Nacelles are assumed to swallow air supersonically. That is, the duct 
is operating at a mass flow ratio of unity. Consistent with this 
assumption, the equivalent body cross-sectional area distribution is 
increased by the oblique projected duct capture area at all stations ahead 
of the duct which are intercepted by an oblique plane. 

Blunt base components are extended (maintaining constant cross- 
sectional area) sufficiently far downstream to prevent flow closure around 
the base. 

In addition to a geometric description, a definition of the pressure 
distribution acting on the configuration is required. The vortex panel 
analysis is used for this purpose. The thickness pressures for planar 
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components have tacitly been neglected under the assumption that the 
surfaces are sufficiently thin that the net pressure coefficient is 
representative of pressure acting on the oblique section. 

Estimation of the wave drag based on equation 20 depends on solution of 
integrals of the type 

1 1 

1 “ | J G "(x 1 )G"(x 2 )ln|x 1 - x 2 |dx 1 dx 2 

0 0 


of a numerically given function G(x) . Evaluation of such forms has been 

16 17 

studied by Eminton ' for functions having G'(x) continuous on the interval 
(0,1) and G'(0) - G'(l) = 0. In such situations, G'(x) can be expanded in a 
Fourier sine series. It can then be shown that 

x - 

N-l 

where 


a n = 


7T 

J G' (x)sin(N^)dV> 
0 


Eminton then solved for the value of the Fourier coefficients which result 
in I being a minimum, subject to the condition that the resulting series for 
G(x) be exact for an arbitrarily specified set of points (0,1), e., i«l,n. 

This approach produces the following result 

n n 

*y~ y~ c.c.f. . 

I— Z_ i J ij 
i=l j=l 

where 


■j- [G(l) -G(0) ] 2 


C i - G(e ± ) -G(0) -[G(l)-G(0)]/i. 




P. . 


cos ^(l-2c.) 

l 


e . 
l 


(n+1) 
-1 


f . . - [P. . ] 
ij 1 ij J 


(£ i' € j )2 ln 


2 (l-2e . ) I 

l 

y'TTarry , 

1 < i < n 


jVlL 2 Vi + 2ye lC ,(l-e.)(l- £ ,) ] 


[€ i +€ j- 2 Vj- 2 ^ £ i € j (1 * € i ) < 1 - e j> 1 J 


+ 2 


(£ i +£ j - £ i £ j >^ i £ j < 1 - £ i )< 1 - e j ) 
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The solution of equation 20 for wave drag is accomplished by use of the 
following identities. 


G(e 8) - A (e 9) 

l e i 

C D ( * )S REF " 1(9)/ L ’< 

w 


2n 



n/2 

- 1 / C (9) M 

n -n/2 W 
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DRAG DUE TO LIFT 


In the discussion which follows, certain limiting (zero and one hundred 
percent suction), and attainable edge force polars are defined. They are 
related to one another as indicated in the following sketch. 


If CD100 > CDO, there is a 

slenderness problem in the 

CD100 calculation. In this 

case the CDO calculation 

is considered more reliable 

and the CDL curve is found 

relative to CDO and the suction level. 



The fixed one hundred percent suction drag due to lift (i.e. 


given by equations 20 through 22. Specifically 



is 


C 


D 


100 



M<1 


+ CL 


V 


w 



THICKNESS 


M>1 


The zero suction drag due to lift is calculated by numerically 
integrating the net pressure distribution times the projected area in the 
streamwise direction over each of the planer surfaces. The following block 
integration scheme is used to sum over all quadrilateral panels. 


where 


and 


'P. 

l 


N 


F l< S REF>L C P. A i“i 


'p. 

1 


i-1 

+ a 3CL /dot + 6 dC /d& 


=0 


PL 

i 


a n . + a + a . 6 


'i i 

a 0 ^is due to twist and camber, S is the control surface deflection, and o ^ 1 
for control surface panels and a ** 0 for non control surface panels. F x ** 2 
for symmetric geometries and F x — 1 for asymmetric geometries. 


Edge forces are neglected in this evaluation. 

The drag due to lift for the total configuration is based on linearized 
potential (100 percent leading edge suction) calculations plus corrections 
to account for suction losses and associated edge vortex forces 
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Ns 


100 


V “ C D_ + S REF ) (1 ' K s. )C s. C j As j°D. + C T (S REF ) 1 } C s, A(x/c T } k T D 

i»l J J J i. _i k k 


Nc„ 


Ns 


C + S 
0 


1 \ 

IEF / 


k-1 

Nc„ 


D„ ■ “REF / (1 ' K s. )C s. C i As i n D. + C T^ S RF,F^ } C c A(x/ Crr ),J 


'»/ *, j j B. T v REF' 

i-1 J J J 


k-1 


s. ' ' I k D. 
k k 


and the leading edge and side edge rotation factors, and T , are (see 
derivation in Appendix B) 


a. 


cos a cos A cos 5 + sin a (cos 6 sin S - sin 8 sin A cos 6 ) 


+ A o/ I A o I [ -cos a cos A sin 8 + sin a (cos 8 cos 8 + sin 8 sin A sin 5 ) ] 

where 6 is the slope angle of the camber line perpendicular to the leading 
edge and the sign of coefficient A 0 from equation 13 is used to determined 

the direction of vortex rotation. 


Tp - ± sin a sin 


C |C | (-cos a sin 8 + sin a cos 8 cos 
n 0 n 0 


S ) 


where 8 is the chordwise slope angle of the tip camber line, plus refers to 
the left side and negative to the right side of the configuration and the 

sign of coefficient C from equation 14 is used to determine the direction 

n o 

of vortex rotation. 


An estimate of the average level of leading edge suction for the 
complete configuration is based on the following equation: 


SUCTION 


< C n * )/(C n -C r 


) , where C 


for K 


100 


0 


= 1.0 ; — > L.E. Suction, if any, is totally recovered. 

^ 1-0 ; => L.E. Suction is partly recovered, the 

remainder is converted to vortex lift and 
drag. 
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HYPERSONIC 


High Mach number analysis is based on non- interfering constant pressure 

1 8 

finite element analysis. 

An arbitrary configuration is approximated by a system of plane 
quadrilateral panels as indicated in figure 13. 



Figure 13. Configuration Represented by Surface Quadrilateral Panels. 


The pressure acting on each panel of a vehicle component is evaluated 
by a specified compres s ion- expans ion method selected from the following 
options . 
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Impact Flow 


Shadow Flow 


1 . 

Modified Newtonian 

i. 

Newtonian (C - 0) 
P 

2. 

Modified Newtonian+Prandtl -Meyer 

2. 

Modified Newtonian+Prandtl -Meyer 

3. 

Tangent wedge 

3. 

Prandtl-Meyer from free -stream 

4. 

Tangent -wedge empirical 

4. 

OSU blunt body empirical 

5. 

Tangent-cone empirical 

5. 

Van Dyke Unified 

6. 

OSU blunt body empirical 

6. 

High Mach base pressure 

7. 

Van Dyke Unified 

7. 

Shock - expansion 

8. 

Blunt -body shear force 

8. 

Input pressure coefficient 

9. 

Shock- expansion 

9. 

Free molecular flow 

10. 

Free molecular flow 



11. 

Input pressure coefficient 



12. 

Hankey flat- surface empirical 



13. 

Delta wing empirical 



14. 

Dahlem-Buck empirical 



15. 

Blast wave 



16. 

Modified tangent-cone 




A discussion of the various methods 

is presented in appendix C. 


Specific analysis recommendations are provided by the program on a component 
by component basis. 

In each method, the only geometric parameter required for determining 
panel pressure is the impact angle, 6 , that the quadrilateral makes with 
the free -stream flow or the change in angle of a panel from a previous point 
where 

S - n/2-0 

cos 6 -(n*V)/(|n| |V|) 

and 

n-ni+ni + n k 
x y^ z 

V = V - Q X r 


H ► -+ — ► -> 

- (V^cos a cos f))i - (V^sin /})j + (V^sin a cos /3)k 
H = pi - qj - rk 

r - (x-x CG )i + (y-y CG )j + (z-z CG )k 

Panel switching between impact or shadow conditions is based on 6 > 0 in the 
former case and S < 0 in the latter. 
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AERODYNAMIC CHARACTERISTICS 


The pressure on each panel is calculated independent of all other 
panels (except the shock-expansion method). If the vehicle is rotating, the 
local pressure coefficient must be corrected to free-stream conditions. 
That is 


s - s, j v i / v . 

local 

Vehicle component forces, which are in the body axis system, are obtained by 
summing panel forces 

AC x “ S REF ^ C p n x A 

AC y“ S REf} C pV 

AC z " S REF ^ C p n z A 

AC i " (bS REF 51 { } C P (z ' Z CC )n y A + ) C P (y ' y CG )n z A } 

AC m “ -^W 1 { } C P (x - X CC )n z A + ) C P (z - Z CC )n x A } 

AC n “ (bS REF^ 1 { } C p (X - X CC )n y A " } C P (y - y CG )n X A } 

where } 

A = panel area 

x,y,z ** coordinates of panel centroid 

igurat ion buildup and total vehicle coefficients are obtained by 
appropriate summation of component contributions. 

The conversion from the body axis system to the wind axis system for 
the lift and drag coefficients is based on the standard trigonometric 
relations. 


c n = C cosa cos/9 - C sinyS + C sina cos 8 
u a y z 

C, — -C sina + C cosa 

^ X 2 

The vehicle static stability derivatives, which are in the body axis system, 
are calculated by the method of small perturbations. Since the basic force 
and moment characteristics are non-linear, these parameters vary with 
attitude angle 
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(Axial) 


C 


x 

a 


C 


z 

a 


C 


m 

a 



C 


C 


n. 


[(C) - (C ) ] / Aar 

a+Aa a 

[(C) - (C ) ] / Aa 

a+Aa a 

I< C „> ' < C m> 

a+Aa a 

[(C) - (C ) ]/Lp 

/3+A/S p 

[(C) - (C ) ] / A/) 

/9+A/3 n p 

[(C) - (C ) ]/A0 

p 


(Normal) 

(Pitch) 

(Side) 

(Yaw) 

(Roll) 


The damping derivatives due to vehicle rotation rate are obtained in a 
similar manner 


V “ l [(C rn ) - (C m ) ]/Aq|> /[(c)/(2V)] 
q ^ q+Aq q ^ 


etc . 


Similarily the control surface derivatives are 



[(C T ) 

(C T ) ] / AS 

S+AS 

L s 

[(C) 

(C ) ] / AS 

S+AS 

m 5 

[(C) 

(C.) ]/& s 

S+AS 

* s 

[(C) 

(C Y ) ] / AS 

S+AS 

* 5 


etc . 

It is the last term in the numerator of these definitions that are being 
calculated and printed in the program output. 
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CONCLUSIONS 


An aerodynamic configuration evaluation program has been developed and 
implemented on a time sharing system with an interactive graphics terminal 
to maximize responsiveness to the preliminary analysis problem. 

The solution is based on potential theory with edge considerations at 
subsonic/supersonic speeds and impact type finite element analysis at 
hypersonic conditions. Three-dimensional configurations having multiple 
non-planar surfaces of arbitrary planform and bodies of non- circular contour 
may be analyzed. Static, rotary, and control longitudinal and lateral- 
directional characteristics may be generated. 

IBM 3081 computation time of less than one minute of CPU/Mach number at 
subsonic, supersonic or hypersonic conditions for a typical simulation 
indicates that the program provides an efficient analysis for systematically 
performing various aerodynamic configuration tradeoff and evaluation 
studies. PRIME 850 and VAX 11/780 computation times are approximately 
fifteen times longer. 
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APPENDIX A 


SUBSONIC/SUPERSONIC FINITE ELEMENT DERIVATIONS 
VELOCITY PERTURBATION POTENTIAL 


The velocity potential for a point source can be used to obtain expressions 
for the velocity potential induced by source and vortex finite elements. 
Integrations are carried out in the z=0 plane and coordinate transformations 
are used to obtain expressions for constant source and doublet panels having 
arbitrary orientation. Consider a surface S having a unit normal, 


n 



n e + n e 
y y z z 


where ( e^, e^, e are the unit vectors in the (x,y,z) coordinate system. 
The velocity potential for a source located at the point (x 0 ,y 0 ,z 0 ) on S is 
given by the expression, 


tf s (x, y,z) 
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2 2 8 d 8 

where □ - (1-M^) — j + — j + — - 

dx dy dz 


The velocity potentials induced by a distribution of sources on the surface 
S a 1 a derived more easily if we transform variables to a coordinate system 

(x,y,z) which has the source distribution on the z - 0 plane. This 
transformation should also preserve the governing differential equation. 
First rotate the coordinate system by <j> to eliminate the y component of the 
normal. 
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In the resulting coordinate system the normal will have components 

( -sin a , 0, cos a ) 

where n = - sin a 

x 


n 


n 


cos a sin <j> 
cos a cos <f> 


2 

cos a 


2 2 

n + n 

y z 


78 



There finally results the following change of variables 
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written 

in terms of the panel normal, 
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If the points (x,y,z), and (x 0 ,y 0 ,z 0 ) both lie in the plane S, then a vector 
joining these points must be perpendicular to n, 

(x-x 0 ) n x + (y-y 0 ) n^ + (z-z 0 ) n z - 0 

A 

and therefore the points lie in the plane, z = 0. This transformation 
preserves the governing Prandtl-Glauert differential equation, since we can 
write , 


2 2 2 2 2 2 
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a-V — + — + — - (1-m j — + — + — 

dx dy dz 3x dy dz 

and if again the point (x 0 ,y 0 ,z 0 ) lies on S, for any (x,y,z), 

A 2 2 A 2 A 2 2 2 2 2 
X + ^ ( y + Z ) - (x-x 0 ) + 0 [ (y-y 0 ) + (z-Z 0 ) ] 
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The velocity potential for an area in the z - 0 plane having constant source 
density is obtained by integrating the influence of infinitesimal source 
elements over the area. Dropping the A and using transformed coordinates we 
have , 
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and since □ <j> 
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- 0 
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A doublet at x 0 , y 0 , 

, z 0 - 

is 
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derivative 

of a 

point source: 
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Integrating from x 0 =£ 0 to infinity yields the potential for a line 
doublet or elementary horseshoe vortex. 
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And an area of constant vortex strength is obtained by integrating this 
expression over the panel area: 
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The solution of these integrals is performed in the following sections. All 
integrals may be checked using tables 1 and 2 at the end of this Appendix. 

The velocity expressions may be obtained by differentiating the 
velocity potentials using table 1. 
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SOURCE PANELS 


First the integration is performed over the panels in the x 0 dire< 
as shown in figure 1. 

x r 
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(x - x x ) - T 21 (y - yi ) 


< X A - x 3 ) - T 4 3 (y B -y 3 ) - o 

Figure 1. Integration Over Panels in x 0 Direction. 


To integrate with respect to y 0 a change of variables is introduced: 


when 
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which is independent of y 0 . 
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Therefore using T b « T + fi , and integrating with £ constant 
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This integration may be checked using table 2. 


Each of the four integration limits corresponds to a corner of the 
quadrilateral. Placing the origin of the x 0 , y 0 coordinate system at one 

2 22 

corner, and setting B - T + , the contribution to $ becomes: 
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and combining each of the four corners: 


BR-(Tx+£ y) 


} 


S>(x,y,z,T) - <f>(,x-x lt 

y-yi. 

2 » ^2 1 ) “ ^ ( X “X 2 , 

y-y 2 . 

2 > *^2 1 ) 

- <Mx-x 3 , 

y-y 3 . 

2 » ^43) *<x-x 4> 

y-y 4 . 

2 » ^43) 


VORTEX PANELS 


Analogous to the source panels the integration is first performed in 


the x 0 direction. 
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changing variables and integrating with respect to r\ 
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therefore for one corner or integration limit 
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Velocity Component Transformations 


The velocity expressions may be obtained by differentiating the 
velocity potentials. The results of this are given on page 33 . Since all 
integrations were done in the A coordinate system, we must consider the 
variable transformation to obtain the actual perturbation velocities. 
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To establish that these are the correct perturbation velocities the 
following criteria must be met: 


1. Laplace's equation must be satisfied 
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8 $ + $ + $ -0 
xx yy zz 


or the equivalent u = v 
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u — w 
z z 
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The correct discontinuity or jump in the perturbation velocity must 
occur at the surface of the quadrilateral panel area. For the source 
panel the jump occurs in the normal or w velocity and on the vortex 
panel there must be a jump of constant magnitude in the u perturbation 
velocity over the panel area. The perturbation velocities should be 
continuous elsewhere, except on the trailing vortex sheet of the vortex 
panel . 

The perturbation velocities must go to zero as upstream infinity is 
approached. 

For the vortex panel the trailing vorticity must extend straight back 
to downstream infinity. This means that any discontinuity in the v 
velocity must be zero outside the spanwise boundaries of the panel and 
must be zero upstream of the panel. 


The first criteria can be established by using the derivatives given in 
table 1. 

The second criteria can be established by noting that all terms except 


tan ^ 2 and tan x 

xy-T(y +z ) z 

are continuous at z - 0. Consider these terms keeping in mind that the 
contributions from all four corners must be included. 


If we let 


£ = (x- Xl ) - T(y-y^ 0 = (x-x ) 5 T(y-y ) 2 

2 2 2 2 2 
R i “ (x * x i> + P Ky-yp + z ] 

and use 

- 1 _i _ i a 4- R 

tan A + tan B = tan ~ 

1 - AB 

then the contributions from both corners on the leading edge can be combined 
as follows . 
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If we define 
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The discontinuity, or jump, in f(z) at z - 0 becomes, 
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Af(z) - lim [f(z)J - lim [f(z)] 
z-+0 z->0 



Therefore when a similar procedure is carried out for the trailing edge 
of a source panel and we subtract the results, we obtain the following jump 
in the w perturbation velocity. 
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For the vortex panel (subsonic) we have an additional term. Considering 
both additional terms from the leading edge corners: 
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after summing all four panels corners (both edges), we obtain the following 
for Au and Av 



|<- region of -+| 

| trailing | 

| vorticity | 

i i 

y=yi=y 3 


To verify the third criteria we must show that all of the functions 
approach zero when all four corners are considered as x -* 
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Therefore considering both corners on the leading edge of the panel 
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and therefore this limit is also zero when both corners of the leading or 
trailing edges are considered. Since all terms are accounted for, the 
perturbation velocities are zero far upstream. 


22 222 222 

Since B R - (Tx + P y) = p [(x-Ty) + B z ] 


2 2 2222 
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there is an apparent singularity along the line 


(x-Ty) * 0 , z =» 0 


However this singularity may be removed by combining the contributions from 
both corners of the leading or trailing edges of the panel. Along either of 
these edges the values of 


(x-x t ) - T(y-y i ) and z 


are the same for each of the panel corners. 
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It can be seen from the above diagram that (Tx 4* fi y) will have the same 
sign on a point (x,y,0) which lies outside the spanwise boundaries of the 
quadrilateral. Therefore outside the spanwise boundaries the term 


2 22 

log [(x-Ty) + B z ] 


can be canceled by combining both corners, and the resulting term 


1 , BR x ±[T(x- Xl )+/? (y-y x ) ] 

± - j lo S 2 

b br 2 ±[T(x-x 2 )+^ (y-y 2 ) ] 


will not be singular if the correct + or - sign is chosen. Within the 
spanwise boundary an actual singularity occurs on the panel edge. 


i 

The term j 
written 


log 
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also has a possible singularity. 


This term can be 
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For the source panel the singularity may be removed for points along 
2 2 

y + z =0 which are outside of the panel boundaries. 
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If (x-x t ) and (x-x 3 ) have the same sign the combination of the two terms 
gives 

! R^X-Xj) ! R 3 + (x-X 3 ) ! R^X-Xj) 

2 log - 2 log “ ± 2 lo S 

R^x-Xj) R 3 -(x-x 3 ) R 3 ±(x-x 3 ) 

where the correct sign is chosen to remove the singularity. On the panel 
edge the singularity is real and cannot be removed. 
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For a vortex panel the terms (subsonic) 
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Both have real singularities for x > 0 (downstream) and removable 
singularities for x < 0 (upstream) . The real singularities occur on the 
panel edges and on the edge of the trailing vortex sheet. 
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and combining the contributions from leading and trailing edges, 
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VELOCITY FLUX FROM AN INCLINED (BODY) SOURCE PANEL 


The perturbation velocity normal to a panel surface is given by the 
expression 
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Since across the panel surface A u = 0, the rate of outflow from the panel 
surface is given by, 
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and since Aw = a 
across the panel surface 
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where 4- and - signify the upper and lower surfaces. 
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SUPERSONIC VELOCITIES - SPECIAL CONSIDERATIONS 


The velocity perturbation influence equations for supersonic flows are 
treated by taking only the real parts of the expressions. This means that 

2 222 1/2 

R“[x+/3(y+z)] is set equal to zero for points which lie outside 

the downstream Mach cone from any given corner. Therefore, R and 
i R+x 

T l°g are zero f° r points which lie outside the downstream Mach cone. 

2 2 2 

For B — (T + P ) > 0, there are no problems using this method. 
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[ t 2 + p ] 


1/2 


i B where i - / - 1 


i B 2 


log 


(Tx+P y)+iBR 


(Tx+/3 y)-iBr 
and combining two corners, 


1,-1 B R 
tan ~ 

(Tx+0 y) 


1 j [T (x-x x )+p (y-y x ) ]+BRj 

- 2 log ^ 

B [T(x-Xj )+p (y-y^j-BR, 


1 x [T(x-x 2 )+^ (y-y 2 ) ]+br 2 
- 2 log i 

b [T(x-x 2 )+£ (y-y 2 ) ] -br 2 


_L_ „ -1 

tan 


r 2 2 

B | [T(x-x 2 )+/? (y-y 2 ) ] Rj - [T(x-Xj )+p (y-y x ) ] R 2 


} 


[T(x-x 2 )+/3 2 (y-y 2 )] [T(x-x 1 )+^ 2 (y-y 1 )] - B 2 R X R 2 


If z - 0 and either R x or R 2 is zero and we allow the other to approach 
zero, the value of F 2 becomes 


f 2 - 


n_ 

B 

0 


[T(x-Xj ) - p (y-y x ) ] [T(x-x 2 ) - p (y-y 2 ) ] < 0 
[T(x-Xj ) - ^ 2 (y- yi )][T(x-x 2 ) - p (y-y x ) ] > 0 


Therefore if R x and R 2 are zero but we are inside the envelope of Mach cones 
from the leading edge (see figure 2) , the value of F 2 is set equal to 


f 2 


n 


B 


if 


' [T(x-x 1 ) - p 2 ( y- yi )][T(x-x 2 ) - p 2 ( y-y 2 )] < 0 

2 2 
- R x < 0 R 2 > 0 

2 2 2 2 

l (x-Ty) > (P - T )z 
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T = tan A 


[(x-x x ) - TCy-yj)] - [(x-x 2 ) - T(y-y 2 )] > (/3 -T )z 


2 2 2 2 2 2 
P - T >0 (x-Xj) < p [(y-yj) + z ] 

(x-x 2 ) < p [ (y-y 2 ) + z ] 


Z 

* 



Figure 2b Supersonic Leading Edge Mach Cone Envelope 



y 



Figure 2b Supersonic Leading Edge Mach Cone Envelope 


The intersection of the lines determind by 

2 2 2 2 
(x-Ty) - (/3 -T )z 


and 


2 2 2 2 
X - p (y + z ) 


occurs on the line 


[ y - ax 1 

\ z — bx f 


therefore 


2 2 2 
1 - 0 (a + b ) 


2 2 2 2 2 
1 - 2aT+T a - (p -T )b 

2 22 2222 2222 
p ( 1 - 2aT+T a ) = (1-0 a ) (0 -T ) - p (0 -T )b 


4 2 


0 a - 2^9 Ta 

+ T 

- 0 => 

a = ±T 


4 2 

2 2 

2 


P b 

- p - T => 

p b = ± 

therefore the 

line 

is determined by 

2 


2 2 2 

2 2 

Tx - B y 

- 0 ; 

; 2 - W - 

T )x 

or 

2 


2 2 

2 2 

Tx - B y 

- 0 ; 

; (x-Ty) = (y9 - 

T )z 


2 2 
P - T 


1 1/2 
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lim 

T->j9 


2 2 

As T -* j3 (sonic leading edge) the value of (T -/?)-+ 0. In this case 


(Ri-R 2 ) 

t [ (x-Xj ) -T(y-y x ) ] 


[(x- Xl )-T(y- yi )] > 0 
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SONIC EDGES 


2 2 2 
As B - T + p 

the function, 


0, numerical difficulties arise in the evaluation of 


1 ! BR+ ( Tx+yS y) 

- J log — 

B BR-(Tx+£ y) 


1 ! BR+[T(x-Ty)+B y] 

- 2 l°g ~ 

B BR- [T(x-Ty)+B y] 


However, for small values of B , this function can be easily evaluated 
numerically by using a few terms of a series expansion. To generate the 
series, first we set 


2 

a - T(x-Ty)+B y 


b 

8 


2 2 2 2 2 
- ( T - B ) [ (x-Ty) + B z ] 

2 2 
B R 

b 


and therefore 


2 2 


2 2 2 
a B R 


B R + b = b ( 1 + 6 ) 


b 6(1 + 8) 


therefore 

2 

! BR+[T(x-Ty)+B y] 

2 !°g 2 

BR-[T(x-Ty)+B y] 

6 


U 


0 


dt 

172 


- 1/2 


[ 5(1+5) ] 1/2 [ 1 


log 


(W) 1 /^ * 1/2 


2 

1 

3 • 2 
1 

3*2 

1 


t + 


5 + 


S + 


S + 


dt 


(1+6) 

1/2. 

s 1 / 2 

2 

J 

0 

1.3 

1 

2 

1 . 3.5 

1 

2 

2 

2! 

t 

3 

2 

3! 

1-3 

1 

2 

£ 

1 . 3.5 

1 

2 

5*2 

2! 

3 

7*2 

3! 

1-3 

1 

2 

5 

1*3.5 

1 

2 

5*2 

2! 

3 

7*2 

3! 

1*3 

1 

2 

£ 

1.3.5 

1 

2 

2 

2! 

3 

2 

3! 


( t(l+t) ] 1/2 


3 

t + 


s + 


3 

s + 


3 

s + 
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a B R 
b 



2 

8 + 
3 


3 

2 


3*5 


2 

s 


4 

2 


5*7 


3 

5 + 


7 

2 


5 * 7*9 


4 

s 


8 

2 


7 * 9*11 


5 

8 


+ 


10 

2 


3 * 7 * 11*13 


6 

8 


1 1 

2 


5 * 9 * 11*13 


7 

6 


+ 


1 5 

2 


5 * 9 * 11 * 13*17 


8 

8 


16 

2 


5 * 11 * 13 * 17*19 


9 

8 


2 


1 8 


+ 

3 * 7 * 11 * 13 * 17*19 


1 o 

8 


1 9 

2 


3 * 7 * 13 * 17 * 19*23 


i i 

8 


+ 


22 

2 


7 * 13 * 17 * 19 * 23*25 


1 2 

6 


2 3 

2 


7 * 17 * 19 * 23 * 25*27 


1 3 

6 


2 

This series converges rapidly for small values of 8, or small values of B . 
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TABLE 1 


TABLE OF DERIVATIVES 


2 2 2 2 2 2 2 2 
y + z , R = x + £ r , 0 - 1 - M 


R f" R 

3x 


d_ - ^ 

» 106 


R+x 


R-x 


1 

R 


g_ 1 i BR+(Tx+/? y) 

9x ~ 2 1°6 2 

B BR- (Tx+^9 y) 

a - 1 z R 

- tan — j 

xy-Tr 


xy-Tr 


2 2 2 

R [(x-Ty) + B z ] 

1 z (Tx+/9 2 y) 

_ 2 2 2 
R [ (x-Ty) + B z ] 


R f- R 

ay 


2 

P y 


a_ i 
ay 


R+x 


1 log 


R-x 


1 x BR+(Tx+/3 y) 

- - 7 log — 

y B BR- (Tx+/9 y) 

2_ -» z R 

3y tan 2 

J xy-Tr 


1 xy 

“ ~~2 
R r 

2 2 

1 x(x-Ty) + p z 

~~ 2 2 2 

R [(x-Ty) + B z ] 

1 Tz (Tx+/? 2 y) 

— 2 2 2 

R [(x-Ty) + B z ] 


r f- r 

3z 


2 

P Z 


d_ i , 

az 2 lo s 


R+x 


R-x 


1 xz 
_ — 

R r 


q_ 1 i BR+(Tx+/5 y) 

dz ~ 2 2 

B BR- (Tx+/9 y) 

.i z R 

3z tan — —T 

xy-Tr 


1 z (Tx+/) y) 

R [ (x-Ty) 2 + B 2 z 2 ] 

1 (x-Ty)(Tx+/y) 

_ 2 2 2 
R [(x-Ty) + B z ] 


2 2 2 
B - T + p 


1 xz 

_ ~2 
R r 


1 xy 
— ~ 

R r 
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TABLE 2 


TABLE OF DERIVATIVES 


2 2 
v + r 


2 2 2 2 2 
(£+>?) + (b -1 ) P , T b 


R f? R 

d_ i 


R+ (Z+n) 


a? 


2 lo e 


R- <£+>?) 


* i , 

a£ 2 log 


S_ -i 

tan 

a? 


bR+(£+b n) 
bR-(£+b%) 

r R 




1 
R 
1 


2 

r -£«? 


2 2 2 

R [£ +b f ] 

1 f (*+b%) 

~~ 2 2 2 
R [£ +b f ] 


R f— R 

drj 


(C+b f?) 


a_ i 

dn 


2 lo g 


R+(fn?) 

R- (£+»?) 

1 i bR+(£+b%) 

3n 2 ^ og 2 

1 b bR-(£+b i,) 


- 1 

a tan 

dr) 


r R 


e«i-f 


i (r -*i7> 
2 2 ~ 
R («7 +f ) 


1 

R 

i r (€+•?) 
~~ 2 2~ 
r (i? +r ) 


e If r 


(b -1) f 


1 R+(^) 


If * iog 


R- (£+»?> 


1 f (?+i7) 
2 5“ 

R (17 +f ) 


l 1 bR+(e+b n) 

5f 2 1°8 2 

11 b bR-(£+b t)) 


1 f (t+b n) 

~ 2 2 2 

R t€ +b r ] 


d_ -i ^ R 

a tan 


i e ce+b ,) 

~ 2 2 2~ 
r te +b f ] 


2 2 
T + £ 


1 *? (fH?) 
2 2 ~ 
R (17 +r ) 
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APPENDIX B 


SURFACE EDGE FORCES 


LEADING EDGE POTENTIAL SUCTION 

In the limit as the wing thickness goes to zero, the increasingly 
reduced pressure acting over a decreasing area results in a limiting suction 
force at the leading edge. If we consider the leading edge region, (figure 
1) the force on the airfoil may be obtained by integrating over a control 
surface in the flow, 


U 

oo 


X 


S 



Figure 1. Leading Edge Suction Region. 


where S is a control surface into which the leading edge penetrates and F is 
the force on the area enclosed by S . In two dimensions the surface integral 
becomes a line integral and since for incompressible, irrotational flow 


P - P 

co 


] 2 2 
— p (u + V ) 


F 

X 


H 


n dS — dy e - dx e 
x y 


l 2 2 

[ P - ~T~ p (u + v ) ] dy + pu [u dy - v dx] 


} 


--H dy + 

C 


1 f 2 2 

p [2uv dx + (v - u ) dy] 

C 


where C is the contour around the leading edge of the airfoil and F is the 
force per unit of leading edge length. 

As the wing thickness approaches zero, the wing becomes a line segment 
(figure 2) and the flow in the leading edge region is identical to the flow 
around a 180 degree corner. Incompressibly , it is described by 
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u = — 1 — cos(0/2) + U cos 0 

r /- oo 

y r 

u = - ~ sin(0/2) - U sin 6 

t ij— 00 

vr 




Figure 2. Wing Represented By Line Segment. 


where (r,0) is a coordinate system centered at the leading edge, and 

u = u cos 0 - u. sin 6 
r 9 


v = u sin 9 + u„ cos 9 
r 0 


dy — R cos(0) d0 
dx = - R sin(0) d0 
and C is the circle r = R. 


Therefore, since 

| P^dy = 0 


/— A * 

7R u = a [ cos(^) cos 9 + sin(^0) sin 9] - a 


7r v = a [ cos(^0) sin 0 + sin(^0) cos 9] - a 


2tt 

2 r r i i i 21 

F^ = pa - cos(*j0) sin(^0) sin 0 + “ — [ sin (^0) - 

0 


2n 

X 2 r 2 2 2 

~ 2 ~P a [ sin ® + cos 0 ] d0 = - 7r pa (1) 

0 


COS (T0 ) 


sin(T0 ) 


2 1 1 
cos (^0) ] cos 9 Y 6.9 
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To relate -F^ , the leading edge suction force, to the pressure distribution 
near the leading edge, the AC across the line segment must be evaluated. 


On 

the 

top 

II 

o 

a 

u - 

Jx 

+ U 

oo 

1 

O 

On 

the 

bottom 

6 *• 2n 

a 

Jx 

+ u 

oo 

, V - 0 


AP 



U 

oo 


4a 

Jx 


4a 1 

and if c is the chord length AC - — ; £ = x/c 

P U Jc j£ 

oo 


c 


t 


A THRUST 
c Az q 

oo 





2 

4a ] 


/c U 

00 


or AC 

P 




A - 


4a 

U JZ 

oo 


These expressions relate the leading edge thrust coefficient to the net 
distribution , AC^, at the leading edge. 

In general we can write 


Cp(^) = A 0 cot(^) + ^ sin(n^) 


n=l 


where 


cot = 


x 

c 


COS(T^) 


sin 


—(1-cos <j>) 

1/2 


2 1 

sin 


1-sin ( 




sin (T<t>) 


[^] 


1/2 


- 1/2 

A 0 is the coefficient of the £ ' term and therefore determines the leading 

edge suction force since only the term which is infinite at the leading edge 
contributes to the suction. 


C 


t 



M = 0 
00 
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For linearized compressible flow the following Mach number correction must 
be applied 

c t - i" * A ° f - 1 - m 1 

To derive the expression for a swept wing, an infinitely skewed wing (figure 
3) is considered 



/ 


Figure 3. Infinitely Skewed Wing Representation. 

Let the subscript or superscript 0 denote the variables normal to the 
leading edge . Then 



the ratio of thrust per unit length is identical in either system 



c C q 

t M co 


At, 

Ay 0 


At 

Ay 


ACp and C^_ in the freestream coordinate system are based on freestream 
dynamic pressure q . Thus 


0 2 
q « q cos $ 

^QO ^OO 



0 2 
AC cos 6 
P 
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0 


C cos 


therefore 


A - A cos 8 
n n 


2 2 2 22 2 2 2 
and /S 0 - 1-M - cos 8 [1/cos 6 -M ] - cos 0 [tan 0 +(1-M ) 


2 2 2 
cos 6 [tan 6 + p ] 


therefore combining terms 


c t (y) 


2 I® C 0 


c o A y q Q 


— b k — — 

8 Po 0 q c 


2 2 

, 2 2 1/0 A 0 COS 8 COS 
cos 6 [ tan 8 + fi ] 4 

8 cos 0 


22 1/2 2 
-f“ [ tan 6 + p ] v A 0 


when AC^ is given by 


AC - A 0 cot(T<£) + 
P 




sin n <f> 


SIDE EDGE POTENTIAL SUCTION 


The method used to compute the suction force at surface tips is similar 
to that used for the leading edge. Since the flow is irrotational 


Av(x,y) 


2 ir AC P (x * y) 


introduce a change of coordinates 

let £ be the fraction of chord 


T be the slope of a constant £ line T - T(£,rj) 
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then integrating 


Av(C^) - — T(£ , »?) AC (i,r,) - 


^ l { [ “ ], + c(,) [fc], 


AC p (|,r,) d£ 


( 2 ) 


Near the tip, we assume a net pressure coefficient of the form 


ACpCC,*?) - 


_1 

r) 

max 



2 

max 



c 

■ ■ av g 
c 


S f(0, 

iN 0 


where 

c 

C N ~c 

avg 

Differentiating 


_1 

r\ 

max 



(V 


max 




1 

J f(Od£ - 1 

0 



ACp (£,»?) 



2 n 

max 



(*7 


max 


2 1-1/2 


avg 


C N f( *> 
LN o 
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Then as r/ -> n , equation 2 gives (keeping only the largest term) 


[ rt ( rj - Y) ) 
max 'max 


- 1/2 


Av(«f,r/) - f t? ( r? - q ) 1 

2 [_ roax max 1 ' J 


- 1/2 


c C 
avg 


N 0 l 


f(x) dx 


U & v (£ >*l) - 2 a ( £ ) rj ( r 7 - >7 ) 

00 [ niax max 7 


-1/2 


-> a(0 - -J- U »?‘ 1/2 c C KT 

o 00 max avg N 0 


| f(x) d X 
0 


Using the expression derived for flow around a corner (equation 1) in 
conjunction with this relation, the suction force at the tip is given by 


AF 


c s (0 


C T AX %> 


2k 1 2 7T avg 2 


o P a 
C T q ® 2 " 


32 


r C 


°T^max 


J f(x) dx 

0 


F 

n 

7T 

1 

C ave _ 2 r 


f f „ , „ 1 

2 ~ 

C T q » 

32 

Vmax N o ] 

4 

| J f(x) dx 

L 0 


dt 


s 

ref ^oo 


c c 
7T avg T 


32 S -n °N 0 
ref max 0 


r r I s 

j J f(x) dx - 

^ 0 


d£ 


where c T is the chord dimension at the tip. 
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EDGE FORCE AND MOMENT INCREMENTS 


To account for edge vortex effects, the linearized forces and moments 
are corrected to reflect losses in suction and the associated formation of 
vortex forces for leading and side edges. The corrections are applied to 
the standard lift, side force and drag coefficients. The corresponding 
increments in the total moment coefficients are calculated by applying the 
above force increments at the appropriate x,y,z coordinates for the leading 
edge stations and center of pressure for the side edges. 

For leading edge force calculations, the lost suction force for each 
span station is given by 


C c As' (1 - K ) 
s s 

where C g is the coefficient of leading edge suction, c is the local chord, 
As' is the local span station width and K g is the leading edge suction 
recovery factor. (K^— 1 - full suction - no vortex) This force is 
subtracted from the direction normal to the section leading edge and re- 
entered as a force component rotated ± 90° about the leading edge. The sign 
of the rotation is determined by the sign of the coefficient A 0 in the 

equation for leading edge suction. 

The change in the total lift, side force and drag is calculated for 
each span station and is written as a function of four coordinate system 
rotations whose rotation angles are known from the leading edge geometry. 

The origin of each coordinate system is located at the leading edge of the 
section camber line. 

The first transformation involves the rotation of the system 
(x 4 ,y 4 , z 4 ) , whose x-axis is tangent to the local normal camber line, to the 

system (x 3 ,y 3 ,z 3 ), whose x-axis it tangent to the corresponding chord plane 

as indicated in figure 4: 


( 



Figure 4. Axis of Rotation for First Transformation in Leading Edge Region. 
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where 5 — 
(dz/dx) c 


-tan' -I [ (dz/dx) + (dz/dx) + (dz/dx) )] / cos A 
l c 6 F 

is streamwise slope due to camber 


} 


(dz/dx) is streamwise slope due to twist 

(dz/dx) is streamwise slope due to flap deflection 

5 F 

A is the local leading edge sweep angle. 

The sweep term converts the total streamwise slope to a slope measured in 
the direction normal to the leading edge. 


The two coordinates systems are related by the following 


matrix: 



cosfi 

0 

-sin5 


o 

0 

1 

0 

< 

c 

y 4 • 

sin£ 

0 

cos5 


c 

Z 4 

L - 


transformation 


The second transformation involves the rotation of the system 
(x 3( y 3i Z 3 ), whose y-axis is tangent to the leading edge, to the system 

(x 2 ,y 2 .z 2 ) • whose y-axis is normal to the configuration center line and in 


the plane of the surface (figure 5) . 



y 2 


s 


Figure 5. Axis of Rotation for Second Transformation in Leading Edge Region. 
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The two coordinates systems are related by the following transformation 
matrix: 


1 


X 2 


c 


y 2 

K 

c 


Z 2 

t j 



0 


0 



fc 1 


X 3 


c 

A 

y 3 


c 


Z 3 


The third transformation (figure 6) involves the rotation of the system 
( x 2 >y 2 > z 2 )i whose z axis is normal to the local surface plane, to the 

system (x 1 ,y 1 ,z 1 ), whose x, y 2 and z axes are in the body axes direction. 



Figure 6. Axis of Rotation for Third Transformation in Leading Edge Region. 


The rotation is about the x 2 , x x axis and of magnitude 6 , the local dihedral 

angle. The two coordinate systems are related by the following 
transformation matrix: 


X 1 


10 0 


c? 

X 

IO 

Vi 


0 cos# -sin# 

« 

c 

y 2 • 

' z i 


0 sin# cos# 


c 

z 2 


The fourth and final transformation (figure 7) involves the rotation of 
the body axis system (x 1 ,y 1 ,z 1 ) to the wind axes system (D , Y, L) . 



D 

x i 


Figure 7. Axis of Rotation for Fourth Transformation in Leading Edge Region. 
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The rotation is about the (y 1 ,Y)-axis and of magnitude a , the angle of 
attack. The coordinate systems are related by the following transformation 
matrix : 




cosq 0 sina 


fc 1 
X 1 

■ C Y 

► = 

0 10 

< 

c 

y i • 

C L 


-sina 0 cosa 


c 

z i 







The composite transformation between the (x 4 ,y 4 ,z 4 ) coordinate system and 
(D,Y,L) coordinate system can then be expressed as 



where is the rotation matrix obtained from multiplication of the four 
previous specified transformation matrices. 


Expressing C ,C and C , in terms of the leading edge suction 
x 4 y* z 4 

parameters , 

C = C c As' (1-K ) 
x 4 s s 


C - 0 

y 4 


C - A /I A I c c As' (1-K ) 
z 4 o' 1 o 1 s s 


we can now write the change in drag, side force and lift resulting from the 
force rotation at each span station: 


A C D - C s c As' (1-K s ) fl D 

A C - C c As' (1-K ) O. 

Y s si 

A C T * C c As' (1-K ) H 

L s s L 


where 

— [cos a (cos A cos 5) + sin a (-sin 
+A 0 / | A 0 | [ -cos a (cos A sin 8) + sin a 


6 sin A cos 5 + cos 6 sin 8)] 

(sin 6 sin A sin 5 + cos $ cos 6)] 
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- [cos 8 sin A cos S + sin $ sin 6 ] 


+A 0 /|A 0 |[cos 8 sin A sin 6 - sin 6 cos 6 ] 

and 

“ [-sin a (cos A cos S) + cos a (-sin 8 sin A cos S + cos 8 sin 6)] 

+A 0 / | A 0 | [ -sin a(‘Cos A sin S) + cos a (sin 8 sin A sin 6 + cos 8 cos 6)] 

For side edge force calculations, the lost suction force at each chord 
station is given by 

C g c T A(x/c) 

where C is the coefficient of side edge suction, c_ is the tip chord and 
s X 

A(x/c) is the local nondimensional chord increment over which is acting. 
This force is subtracted from the direction normal to the tip chord and re- 
entered as a force component rotated ± 90° about the tip chord. The sign of 

the rotation is determined by the sign of the coefficient C V7 in the 

N 0 

equation for side edge suction. 

In a manner similar to that for the leading edge forces, the change in 
the total lift, side force and drag coefficients is calculated for each 
chord increment and is written as a function of three coordinate system 
rotations whose angles are known from the tip geometry. The origin of each 
coordinate system is located on the chord line at the beginning of each 
chord increment. 

The first transformation (figure 8) involves the rotation of the system 
( x 3 »y 3 » z 3 )» whose X axis is parallel to the local camber line, to the system 

(x 2 ,y 2 ,z 2 )» w h° se axis is tangent to the tip chord. 



Figure 8. Axis of Rotation for First Transformation Along Chord. 

where 8 - tan _1 [ (dz/dx) + (dz/dx) + (dz/dx) f )] 

c € 6 f 
( dz/dx)^ is streamwise slope due to camber 
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and 


(dz/dx)^ is streamwise slope due to twist 

(dz/dx)^ ) is streamwise slope due to flap deflection 
F 

The two coordinate systems are related by the following transformation 
matrix: 

cos<5 0 -sin$ 

0 10 
sin<5 0 cos<S 

The second transformation (figure 9) involves the rotation of the 
system (x 2 ,y 2 ,z 2 ), whose y-axis is normal to the tip chord, to the system 

( x i ,Yi ,Zi) , whose x, y, and z-axes are in the body axes direction. 





Figure 9. Axis of Rotation for Second Transformation Along Chord. 

The rotation is about the (x 2 ,y 1 )-axis and of magnitude #, the local 

dihedral angle. The two coordinate systems are related by the following 
transformation matrix: 

10 0 
0 cos# -sin# 

0 sin# cos# 
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The third and final transformation (figure 10) involves the rotation of 
the body axes system (x 1> y lf z 1 ) to the wind axes system (D,Y,L). 



Figure 10. Axis of Rotation for Third Transformation Along Chord. 

The rotation is about the (y lf Y)-axis and of magnitude a , the angle of 

attack. The two coordinate systems are related by the following 
transformation matrix: 



The transformation between the (x 3 ,y 3J Z 3 ) coordinate system and the 
(D,Y,L) coordinate system can then be expressed as 



where T is the rotation matrix obtained from multiplication of the three 
previously specified transformation matrices. 

Expressing C , C and C in terms of the side edge suction 

x 3 y 3 Z 3 

parameters , 
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C - 0 

x. 


C y, " C s c T A(x/C) 


C z 3 - < V IC No l)C = C T a<X/C) 

we can now write the change in drag, side force and lift resulting from the 
force rotation at each side edge station: 

4C D - C s C T l4<X ' /c>T D 

“ Y - C s c T A(it/c)T Y 
ac l - c s c i 2a < x / c > t l 

where 

T D “ ± sin(a)sin(0) + (C N /|C N | ) [ -cos(a)sin(5) + sin(a)cos(0)cos(6) ] 

T y - ± cos (0) + (Cjj /IC |)[sin(tf)cos(«)] 

T T - ± cos(a)sin(0) + (C / | C | ) [sin(a)sin(S )+ cos(a)cos(0)cos(5) ] 

The minus sign on the first term of each equation is for the right side of 
the configuration and the positive sign is for the left side. These force 
increments are numerically integrated along each tip chord to obtain the 
total change in lift, side force and drag due to side edge force rotation. 
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APPENDIX C 


HYPERSONIC FINITE ELEMENT ANALYSIS 


High Mach number analysis has a number of optional methods for calculating 
the pressure coefficient. In each method the only geometric parameter 
required is the element impact angle, 5, or the change in the angle of an 
element from a previous point. 

The methods to be used in calculating the pressure in impact (5 > 0) and 
shadow (6 < 0) regions may be specified independently. A summary of the 
program pressure options is presented below. 


Impact Flow Shadow Flow 


1 . 

Modified Newtonian 

i. 

Newtonian (C - 0) 

V 

2. 

Modified Newtonian+Prandtl -Meyer 

2. 

r 

Modified Newtonian+Prandtl -Meyer 

3. 

Tangent wedge 

3. 

Prandtl -Meyer from free -stream 

4. 

Tangent-wedge empirical 

4. 

OSU blunt body empirical 

5. 

Tangent-cone empirical 

5. 

Van Dyke Unified 

6. 

OSU blunt body empirical 

6. 

High Mach base pressure 

7. 

Van Dyke Unified 

7. 

Shock- expansion 

8. 

Blunt-body skin friction model 

8. 

Input pressure coefficient 

9. 

Shock - expans ion 

9. 

Free molecular flow 

10. 

Free molecular flow 



11. 

Input pressure coefficient 



12. 

Hankey flat- surface empirical 



13. 

Delta wing empirical 



14. 

Dahlem-Buck empirical 



15. 

Blast wave 



16. 

Modified tangent-cone 




and C are in the stability axis system. Other coefficients are in the 

body reference coordinated system. It should also be noted that side force 
and pitching moment coefficients are invariant in an (a, fi) transformation, 
whereas the yawing and rolling moment coefficents are not invariant. 

A brief review of these methods will be presented in the following text. 


MODIFIED NEWTONIAN 


This method is probably the most widely used of all the hypersonic force 
analysis techniques. The major reason for this is its simplicity. Like all 
the force calculation methods, however, its validity in any particular 
application depends upon the flight condition and the shape of the vehicle 
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or component being considered. Its most general application is for blunt 
shapes at high hypersonic speed. The usual form of the modified Newtonian 
pressure coefficient is 

2 

C - K sin 8 
P 


In true Newtonian flow (M - « , 7 - 1) the parameter K is taken as 2 . In 
the various forms of modified Newtonian theory, K is given values other than 
2 depending on the type of modified Newtonian theory used. K is frequently 
taken as being equal to the stagnation pressure coefficient. In other forms 
it is determined by the following relationship (Reference 19) . 


K 


where 


2 


C 

p 

nose 


/sin 


8 

nose 


C 


P 


nose 


5 

nose 


the exact value of the pressure 
coefficient at the nose or leading 
edge 

impact angle at the nose or leading 
edge 


In other work K is determined purely on an empirical basis. 

K - fn (M, a, shape) 

When modified Newtonian theory is used, the pressure coefficient in shadow 
regions ( 8 < 0 ) is usually set equal to zero. 


MODIFIED NEWTONIAN PLUS PRANDTL - MEYER 


This method, described as the blunt body Newtonian + Prandtl -Meyer 
technique, is based on the analysis presented by Kaufman in Reference 20. 
The flow model used in this method assumes a blunt body with a detached 
shock, followed by an expansion around the body to supersonic conditions. 
This method uses a combination of modified Newtonian and Prandtl -Meyer 
expansion theory. Modified Newtonian theory is used along the body until a 
point is reached where both the pressure and the pressure gradients match 
those that would be calculated by a continuing Prandtl -Meyer expansion. 

The calculation procedure derived for determining the pressure coefficient 
using the blunt body Newtonian + Prandtl -Meyer technique is outlined below. 

1. Calculate free-stream static to stagnation pressure ratio 


r 2 l(7/(7-l ))f 2 a(1/(7-D) 

P = PyP 0 - |2/[(7 + 1) Mjj |[2 7 M oo - (7 - l)]/(7 + l)j 
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2 . 


Assume a starting value of the matching Mach number, M 
assume M =1.35) 

q 


(for 7 - 1.4 


3. Calculate matching point to free-stream static pressure ratio 


r 2 ^(7/(7-!)) 

Q - P q /P 0 - |2/[2 + (7 * DM q ]| 

4. Calculate new free-stream static to stagnation pressure ratio 
P c - Q {l - (7 2 M q Q)/[4(M q - 1)(1 - Q)] 

5. Assume a new matching point Mach number (1.75) and repeat the above 
steps to obtain a second set of data. 

6. With the above two tries use a linear interpolation equation to 
estimate a new matching point Mach number. This process is repeated 
until the solution converges. 

7. Calculate the surface slope at the matching point 
sin(5 q ) - (Q - P)/(l - P) 

8. Use the Prandtl-Meyer expansion equations to find the Mach number on 
the surface element, M 

o 

9. Calculate the surface pressure ratio 

p 5 /p 0 = n c [l + (7 - D/<2) m^]* (7/(7 ' 1)) 

where 

rj is provided as an empirical correction factor 

P is the pressure on the element of interest 
o 

10. Calculate the surface to free-stream pressure ratio 

V p „ - < 1 / p ) <v p °> 

11. Calculate the surface pressure coefficient 

Cp - 2/(7MJ(P 6 /P to - 1) 

8 
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The results of typical calculations using the above procedure are shown in 
Figure 1. Note that the calculations give a positive pressure coefficient 
at a zero impact angle. As pointed out in several references these results 
correlate well with test data for blunt shapes. However, if the surface 
curvature changes gradually to zero slope some distance from the blunt 
stagnation point the pressure calculated by this method will be too high. 
This is caused by characteristics near the nose intersecting the curved 
shock system and being reflected back onto the body. If the zero slope is 
reached near the nose (such as in a hemisphere or a cylinder) this effect 
has not had time to occur. 


TANGENT -WEDGE 

The tangent-wedge and tangent-cone theories are frequently used to calculate 
the pressures on two-dimensional bodies and bodies of revolution, 
respectively. These methods are really empirical in nature since they have 
no firm theoretical basis. They are suggested, however, by the results of 
more exact theories that show that the pressure on a surface in impact flow 
is primarily a function of the local impact angle. In this program the 
tangent-wedge pressures are calculated using the oblique shock relationships 
of NACA TR- 1135 (Reference 21). The basic equation used is the cubic given 
by 

[sin (* s )] 3 + b[sin( 0 s )] 2 + c[sin(^)] + d - 0 
or 

3 2 

R + bR + cR + d = 0 

where 

$ -= shock angle 

6 - wedge angle 

2 2 2 

b - -(M + 2)/M - 7 sin( 6 ) 

c - ( 2 M + 1 )/M + [(7 + 1) /4 + (7 - 1 )/M ] sin(S) 

2 4 

d - cos( 6 )/M 
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The roots of the above cubic equation may be obtained by using the 
trigonometric solution procedure (see Reference 22) as indicated below. 


2 J- p/3 cos (w/3) - b/3 
-2 /-p/3 cos (w/3 + 60 ) - b/3 


-2 7- p/3 cos (w/3 - 60 ) - b/3 

y x - b / 3 

y 2 • b / 3 

y 3 - b/3 


roots of the reduced cubic equation 

-(b 2 /3) + c 

2 (b/3) 3 - bc/3 + d 

cos (w) - -q/ (27- (p/3) ) 

2 

R. - sin(0 ) = roots of the cubic equation 

l s 

The smallest of the three roots corresponds to a decrease in entropy and is 
disregarded. The largest root is also disregarded since it never appears in 
phys ical ac tual i ty . 

For small deflections, the cubic solution becomes very sensitive to 
numerical accuracy; that is, to the number of significant digits carried. 
Since this is dependent on the particular machine employed, an alternate 
procedure is used. 

When the flow deflection angle is equal to or less than 2.0 degrees, the 
following equation is used instead of the above cubic relationships 
(Reference 23): 

sin(0 ) - 1/M 2 + (7 + 1 ) / ( 2 ) S A/M 2 - 1 



where 

y i 

p 

q 
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Once the shock angle is obtained the remaining flow properties may be found 
from the relationship of Reference 21. 


density 


temperature 


PrT P ^ 6M sin(0 s )]/[M sin(0 s ) + 5]j- 


{■ 

{> 


T 0 - T \ [7(M 2 sin(^ g )-l)(M 2 sin(tf s )+5)]/[36M sin(0 s )]j- 


pressure 

coefficient 


{ 


[7M 2 sin(^ s )-l]/6j-/(0.7M 2 ) 


where 


( ), 


conditions behind the shock 


Oblique shock detachment conditions are reached when no solution may be 
found to the above cubic relationships. Under these conditions the program 
uses the Newtonian + Prandtl-Meyer method for continued calculations. 


TANGENT -WEDGE. TANGENT -CONE . AND DELTA WING 
NEWTONIAN EMPIRICAL METHOD 


The tangent -cone and the tangent -wedge (figure 2) Newtonian empirical 
methods used in this program are based on the empirical relationships 
derived below. 



Figure 2. Tangent-Cone and Wedge Notations. 


For wedge flow 

sin(0 g ) = sin(<$ w )/[(l - e)cos($ 8^)] 

where 

c = p/p 2 = (7 - l)/(7 + l){l + 2/1(7 - 1)M m 2 ]} 
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For cone flow (thin shock layer assumption) 

sin(0 ) ~ sin(fi )/[(l - e/2)cos(0 - 8 )] 

S c £» t- 

In the limit as H ■*«,£■= e .. - (7 - l)/(7 + 1) and cos (9 - 8 ) - 1 

lim s 

Therefore 

wedge cone 

sin(0 ) - (7 + l)/2 sin(5 ) sin(0 ) - 2(7 + l )/(7 + 3) sin( 6 c > 

These limiting expressions for 8 may now be compared with the data of TR- 
1135 (Reference 21) at 7 - 7/5 using the following similarity parameters. 

The exact equations contain three variables - 8 8 , and c. Noting that for 

7 ” constant, e = fn (M ) only, the preceding equations may be rewritten in 

the following form: 

wedge cone 

M - M sin(5 )/[(l-«)cos(ff - 8 ) ] M - M sin(5 )/[ (l-e)cos(0 - fi )] 

ns w s w ns w aw 

The parameter ( 8 - 6 ) is approximately constant and independent of M except 

near the shock detachment condition. The equations essentially contain only 

two variables, M and M sin 6 . These are used as coordinates to plot the 
ns 

data for wedge flow shown in Figure 3. A similar plot could be obtained for 
cone flow. From the figure it is seen that the data are nearly normalized 
with the use of these coordinates . 

For rapid calculation we need relationships for M^as a function of M 
sin( 5 ) that satisfy the following requirements: 

1. The effect of shock detachment is neglected 

2. At M sin( 6 ) « 0, M * 1 

ns 

3. The solution asymptotically approaches the M - °o line 

4. Have the correct slope, d[M ]/ d[M sin( 6 )] at M sin(S)— 0 

These conditions lead to equations of the following form 

-(K M'/2) 

wedge M - K M' + e w 

6 ns w 

K - (7 + l)/2 

w 
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cone M 

ns 


K M' 
w 


+ e 


<K c M') 


where 


M' “ M sin(5) 

K 2(7 + l )/(7 + 3) 

c 

These expressions are compared with the data of TR-1135 in Figures 4 and 5. 
The cone data are also shown in Figure 6 with the same scales as in Figure 
3. 

The pressure coefficient may now be obtained by the following relationships 
for a wedge and cone respectively. 

C p - 4/(7 + D(\ s - D/M* 

_ l 

C p - 2sin(S)jl - [<7 - D M n S + 2]/[4( 7 + l)M^ s ]| 

Experimental results have shown the pressure on the centerline of a delta 
wing to be in agreement with two-dimensional theory at small values of the 
similarity parameter (M'<3.0) and with conical flow theory at higher values. 
The previous expressions derived for wedge and cone flows have been combined 
to give these features. The resulting relationships are given below. 


M -KM- - e-'W 2 *' 
ns c 


For 7 - 7/5 


M 

ns 


1.09Msin(S) + e 


(0 ,49Ms in(6) ) 


The similarity parameter relationship for pressure is 

M 2 Cp - (4/(7 + D(\ s - !) 

The shock angle and pressure coefficient calculated from the above equations 
are compared with the experimental results (Reference 28) in Figures 7 and 8 
respectively . 
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Figure 4 


I ~ — I I T I 

1.2 1.6 2.0 2.4 24 

M sin 6 =• M* 


Wedge Flow Shock Angle Empirical Correlation. 
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M sin 8 ~ M' 
00 


Figure 5. Conical Flow Shock Angle Empirical Correlation. 
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Figure 7. Delta Wing Centerline Shock Angle Correlation. 
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Figure 8. Delta Wing Centerline Pressure Coefficient Correlation. 
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OSU BLUNT BODY EMPIRICAL METHOD 


The OSU (Ohio State University) blunt body empirical equation describes the 
pressure distribution about cylinders in supersonic flow. The equation was 
presented in Reference 25 and was stated to match "all the data obtained on 
the cylinders in the present test series with a maximum deviation of 2.5 
percent." The expression used is 

Pj/P - 0.32 + 0.455 cos(0) + 0.195 cos(20) + 0.035 cos(30) - 0.005 cos(40) 
00 

where 

6 — peripheral angle on a cylinder 

o 

(- 0 at the stagnation point) - (90 - 5) 

- surface pressure 

— total pressure rise through normal shock 

00 

The pressure coefficient is calculated from the relationship 

C - [ <P 1 /P t )(P t /PJ - l]/(7M 2 /2) 

* oo CO 

where 

P t /P m - K 7 M 2 /2 + 1 

00 

K * stagnation pressure coefficient — C 

P stag 

p^ - freestream pressure 

7 = ratio of specific heats =1.4 

VAN DYKE UNIFIED METHOD 


This force calculation method is based on the unified supersonic -hypersonic 
small disturbance theory proposed by Van Dyke in Reference 26 as applied to 
basic hypersonic similarity results. The method is useful for thin profile 
shapes and as the name implies extends down to the supersonic speed region. 

The similarity equations that form the basis of this method are derived by 
manipulating the oblique shock relations for hypersonic flow. The basic 
derivations are shown on pages 753 and 754 of Reference 31. The result 
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obtained for a compression surface under the assumption of a small 
deflection angle and large Mach number is (hypersonic similarity equation) . 


C p - S 2 [(7 + l)/2 + y<(7 + l)/2) 2 + 4/H 2 ] 

where H is the hypersonic similarity parameter given by MS . The 
contribution by Van Dyke in Reference 26 suggests that this relationship 
will also be valid in the realm of supersonic linear theory if the 


hypersonic similarity parameter (Jh - 1 )S. This latter parameter is used 
in the calculations for this force option in the arbitrary body program. 

A similar method may also be obtained for a surface in expansion flow with 
no leading edge shock such as on the upper side of an airfoil. The 
resulting equation is 

C p - M2/(7H 2 )]{[(l-(7 - l)H/2)] (27/(7 ' 1) . lj- 


2 

where again H is taken to be (7 m - 1 )S in the unified theory approach. 


SHOCK -EXPANS ION METHOD 


This force calculation method is based on classical shock- expans ion theory 
(see Reference 27). In this method the surface elements are handled in a 
"strip- theory" manner. The characteristics of the first element of each 
longitudinal strip of elements may be calculated by oblique shock theory, by 
conical flow theory, or by a Prandtl -Meyer expansion. Downstream of this 
initial element the forces are calculated by a Prandtl-Meyer expansion. 

By a proper selection of the element orientation the method may be used for 
both wing-like shapes and for more complex body shapes. In this latter case 
the method operates in a hypersonic shock -expans ion theory mode. 


FREE MOLECULAR FLOW METHOD 


At very high altitudes conventional continuum flow theories fail and one 
must begin to consider the general macroscopic mass, force, and energy 
transfer problem at the body surface. This condition occurs when the air is 
sufficiently rarefied so that the mean free path of the molecules is much 
greater than a characteristic body dimension. This condition is known as 
free molecular flow and the method of analysis selected for this program is 
described in Reference 28. This method was also used in Reference 29. The 
equations used were taken from these references and are presented below. 
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Pressure coefficient 


C - s' 2 {[(2-f )/Jv S sin(fi) + f /2 7 t,/T ] e ‘ (S 

p ^ n ' n b « 

+ [<2-f n >(s’.in( t > + 1/2) i f n /2/WvTs sin(5) J [1 . «rf(S sin(«>] 
Shear force coefficient 

r . 2 _ 

C f - [cos(fi)f ]/(/* S)-[e' (S sln < 5 '>+ J w S sin(5)[l + erf (S sin(fi))] 


} 

} 


where 


S - speed ratio = 77/2 


f 

n 


S 


normal momentum accommodation coefficient 
1.0 for Newtonian 

0.0 for completely diffuse reflection) 
impact angle 


o 

- body temperature, K 


- free -stream temperature, K 

x 2 

erf - error function erf (x) = 2/tt J e dx 

o 

f - tangential momentum accommodation coefficient 

0.0 for Newtonian flow 

1.0 for completely diffuse reflection 

The pressure force acts perpendicular to the surface and this direction is 
readily obtained since the element normal has already been determined in the 
geometry subroutines. The shear force acts in the direction of the 
tangential velocity component on the surface and this direction is 
determined by taking successive vector products. The procedure is 
illustrated in figure 9 where the incident velocity vector is defined as 


v - v x i + v y j + v z k 
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and the surface normal as 



First, a surface tangent vector (T) is defined by the cross product of the 
normal and velocity vectors; 


where 


T x 1 + t y 3 + T z k 


N V - N V 
Y Z Z Y 


N V - N V 
Z X X Z 


N V - N V 
X Y Y X 


Then the direction of the shear force (S) is given by the cross product of 
the surface tangent and normal vectors; 


S - s x i + S y j + s z k 

where 

s x - t y n z - T z n y 

S Y ” T Z N X T X N Z 

s z “ T x n y t y n x 

The final components of the shear force in the vehicle axis system are given 
by 
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(SHEAR) (S ) / STOTAL 

A 

(SHEAR) (S ) / STOTAL 
(SHEAR) (S ) / STOTAL 

where 

SHEAR is the shear force as calculated by the free molecular flow 
equations . 

2 2 21/2 

STOTAL - (S v + S v + S_ ) ' 

X Y Z 

In using the free molecular flow method the above analysis must be carried 
out over the entire surface of the shape including the base, shadow regions, 
etc. When the free molecular flow method is selected, it is used for both 
impact and shadow region. 

The plane formed by the velocity vector and the surface normal is referred 
to as the velocity plane (shaded region in the sketch), since both the 
incident and surface velocity are in this plane. This definition is correct 
for two-dimensional flow, however, it is only an approximation to the shear 
direction in the general arbitrary-body case. 


SHEAR X 

SHEARy 

SHEAR„ 


HANKEY flat- SURFACE empirical method 


This method uses an empirical correlation for lower surface pressures on 
blunted flat plates. The method, derived in Reference 30, approximates 
tangent-wedge at low impact angles and approaches Newtonian at high impact 
angles. The pressure coefficient is given by 

2 

0^ - 1.95 sin(S) + 0.21 cos(6) sin(<$) 


DAHLEM-BUCK EMPIRICAL METHOD 


This is an impact method that has been derived such that tangent -cone and 
Newtonian results are approximated, respectively, at low and high values of 
the impact angle. The empirical relationships presented in Reference 31 are 


for 

6 < 

22.5 

for 

S > 

22.5 


C p - + [ sin(46 ) 

C - 2.0 sin(S) 

P 


3/4 


} 


sin(6) 5/4 /[4cos(S)cos(26)] 3 / 4 


135 



BLAST WAVE PRESSURE INCREMENTS 


This method uses conventional blast -wave parameters to calculate the 
overpressure due to bluntness effects. Force contributions determined by 
this procedure must be added to the regular inviscid pressure forces 
(tangent -wedge , tangent-cone, Newtonian, etc.) calculated over the same 
vehicle geometry. The specific blast wave solutions used in the Program 
were derived by Lukasiewicz in Reference 32: 


P/P -AM 


{(C D ) ( 1 /( 1 +j)) /[(X 0 - X)/d]j> 


(2+j )/3 


+ B 


where 


Cp is the nose drag coefficient 

d is the nose diameter or thickness 
X 0 is a coordinate reference point 

and the coefficients A, B are 


Flow 

} 

A 

B 

Two dimensional 

0 

0.121 

0.56 

Axisymmetric 

1 

0.067 

0.44 


MODIFIED TANGENT -CONE METHOD 


This method, originally developed for use on cones with elliptical cross 
sections, modifies the tangent-cone result by an increment representing the 
deviation from an average pressure divided by an average Mach number. More 
specifically, the following equations are used (after Jacobs, Reference 33) 


C 

P 

where 


C - (C - C )/M 
P tc P tc P avg aV S 


tc 


tc 


is the surface pressure coefficient 

is the conventional tangent-cone pressure coefficient 
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is the average pressure coefficient 


C 

P 


avg 


1 


\ k/ } a ' 


A is element area 


M & *- s t ^ ie avera S e Mach number, defined for an equivalent 

cone having pressure coefficient C 

P avg 


HIGH MACH BASE PRESSURES 


For a body in high speed flow it might be expected that any base regions 
would experience total vacuum. That is, 


c p - - V(7My2) 

However, the viscosity of real gases causes some pressure to be felt in base 
region and experimental data have shown this to be roughly 70 % vacuum for 
air. Therefore, the expression 


C 

P 



has been included in the program. 
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